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ZONAL  MODELS  OF  TURBULENCE  AND 
THEIR  APPLICATION  TO  FREE  SHEAR  FLOWS 


Abstract 

The  concept  of  zonal  modeling  Is  tested  by  using  three  homogeneous 
flows.  Turbulence  models  are  constructed  Independently  for  various 
zones.  A  parameter  that  governs  the  readjustment  of  the  flow  from  one 
zone  to  another  is  Introduced.  This  readjustment  parameter  reflects  the 
changes  in  physics  Involved  in  the  transition  process.  The  results  of 
the  zonal  concept  are  very  promising. 

The  idea  is  then  applied  to  the  development  of  Improved  models  for 
free  shear  flows.  A  close  examination  of  the  experimental  data  of  vari¬ 
ous  free  shear  flows  is  first  performed.  The  standard  K-e  model  Is 
adopted  as  a  base  model  since  it  predicts  many  flows  reasonably  well. 
Zonal  models  for  well-defined  zones  are  constructed  Independently  by 
systematically  modifying  the  standard  K-e  model  to  reflect  the  physics 
associated  with  the  zones. 

It  Is  found  that  two  governing  parameters  suffice  to  classify  all 
free  shear  flows  considered,  and  they  are  therefore  used  in  the  present 
model.  >The  first  parameter  represents  the  ratio  of  the  inertial  force 
to  the  driving  force  of  the  shear  layer  or,  alternatively,  the  ratio  of 
eddy  turn-over  time  to  the  characteristic  free  stream  time.  It  deter¬ 
mines  the  relative  Importance  of  the  diffusion  process  in  the  turbulence 
transport  mechanism.  The  larger  this  parameter,  is  the  more  Important 
the  diffusion  process  of  turbulence  becomes.  The  other  parameter  is  a 
measure  of  lateral  vortex  stretching.  Strong  vortex  stretching  reduces 
the  correlation  between  turbulent  shear  stress  and  turbulent  kinetic 
energy  and  hence  reduces  the  spreading  rate  of  the  flow. 

There  are  a  total  of  three  zones  found  in  the  class  of  free  shear 
flows;  they  represent  the  limiting  values  of  the  two  governing  param¬ 
eters.  The  zonal  models  for  these  three  zones  are  built  separately. 
They  are  then  blended  to  allowed  readjustments  between  zones  by  means  of 
readjustment  parameters  that  are  functions  of  the  governing  parameters. 
As  a  result,  a  unified  zonal  model  for  all  free  shear  flows  is  produced. 

This  unified  zonal  model  is  used  to  predict  five  basic  free  shear 
flows.  The  agreement  with  data  Is  excellent.  Comparisons  of  model  per¬ 
formance  are  made  with  the  standard  K-e  model. 
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Chapter  I 


INTRODUCTION 

1.1  General  Introduction  to  Zonal  Modeling 

Methods  of  simulating  turbulent  flows  can  be  classified  according 
to  the  following  scheme  (Kline  et  al.  (1981)): 

1.  Correlations 

2.  Integral  methods 

3.  One-point  closure  methods 

4.  Two-point  closure  methods 

5.  Large-eddy  simulation 

6.  Full  simulation 

As  one  moves  downwards  in  the  list,  each  method  requires  more  computa¬ 
tional  resources  but  less  modeling  than  those  above  it;  consequently, 
one  might  think  it  should  also  be  more  accurate.  However,  the  higher- 
level  simulations  also  require  more  complex  data  inputs  which  are  often 
not  as  accurately  known.  Therefore,  a  simulation  at  a  given  level  is 
not  always  more  accurate  than  simulations  at  lower  levels.  This  remark 
applies  primarily  to  levels  2,  3,  and  perhaps  4.  Also,  owing  to  the 
tremendous  cost  of  simulating  even  the  simplest  flows  by  means  of  high 
level  simulations  (large-eddy  simulation  and  full  simulation)  on  present 
computers,  such  simulations  can  serve  only  as  research  tools  and  not 
engineering  design  tools  at  the  current  time. 

In  the  trade-off  between  accuracy  and  computation  time,  one-point 
closure  methods  seem  to  offer  the  best  compromise  for  high  technology 
applications  at  the  present  time.  The  one-point  closure  category 
includes:  (i)  mixing-length;  (ii)  one-  and  two-equation  models;  and 

(iii)  Reynolds  stress  and  algebraic  stress  models;  these  three  approx¬ 
imations  are  sometimes  referred  to  as  yesterday's,  today's,  and  tomor¬ 
row's  models,  respectively.  The  latter  remark  is  based  on  the  idea  that 
this  ordering  represents  a  ranking  of  the  quality  of  the  models,  i.e., 
that  each  succeeding  model  is  more  exact  than  its  predecessors;  this  may 
not  always  be  the  case  (see,  for  example,  Kline  et  al.  (1981)). 
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On  balancing  all  theae  factors,  we  believe  It  Is  better  to  adopt 
the  pragaatlc  view  that  turbulence  aodels  are  engineering  correlations 
In  a  aore  aophlat lcated  gulae  and  that  the  only  true  teat  of  quality  la 
performance.  Of  the  aodela  mentioned,  Reynolds  streaa  aodels  are  still 
In  the  developaent  stage,  and  two-equation  and  algebraic  stress  aodela 
are  the  aost  popular  onea  at  present  (see,  for  exaaple,  Rodl  (I960)). 
Two-equation  aodels  also  have  laportant  advantages  In  flexibility  and 
range  of  application  coapared  to  one-equation  aodels. 

Two-equation  aodels  consist  of  an  assuaed  functional  fora  for  the 
structure  of  the  Reynolds  stress  tensor  (usually  the  eddy  viscosity 
assuaptlon)  and  two  partial  differential  equations  for  turbulent  kinetic 
energy  and  another  turbulence  quantity  that  provides  a  length  scale. 
The  aost  coaaonly  used  two-equation  aodels  are  K-e,  K-L,  and  K-w  aodels 
(where  K  -  kinetic  energy;  e  •  dissipation  of  turbulence;  L  -  length 
scale;  and  u  -  vortlclty). 

Looking  at  the  results  presented  in  1980-81  AFOSR-HTTM-Stanford 
Conference  on  Coaplex  Turbulent  Flows,  we  find,  for  aany  of  the  flows 
tested,  two-equation  aodels  gave  the  best  predictions.  However,  for 
others,  soae  of  the  slapler  aodels,  including  Integral  and  mixing  length 
aodels,  were  as  good  or  better.  Unfortunately,  which  method  Is  best 
varies  from  case  to  case.  Each  of  the  models  In  the  1980-81  Conference 
required  aodlf lcations  or  special  treatments  for  particular  flows.  No 
single  model  presented  was  accurate  over  the  wide  range  of  cases  used. 
Moreover,  soae  coaputors  reported  success  on  soae  classes  of  flows  fol¬ 
lowed  by  degradation  of  the  results  when  attempts  were  aade  to  extend 
the  range  of  flows  without  change  of  aodel  or  constants. 

Although  having  a  universal  aodel  Is  desirable,  the  preceding  para¬ 
graph  indicates  that  no  such  aodel  exists  at  present.  It  also  appears 
such  a  aodel  would  need  to  be  very  complicated  and  would  probably  be 
difficult  to  apply.  (This  opinion  Is  now  shared  by  W.  C.  Reynolds,  who 
has  been  a  leading  proponent  of  universal  models  for  soae  years,  private 
comunlcatlon. )  It  Is  our  view  that  aodels  should  be  closely  tied  to  the 
structure  of  the  flow  phenomena  they  are  supposed  to  represent.  This 
view  is  shared  by  aany  others  but  has  rarely  been  used  as  the  basis  for 
aodels.  Since  different  parts  of  flow  fields  often  have  quite  different 
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flow  structures.  It  say  be  impossible  to  model  all  flows,  or  even  all  of 
a  single  complex  flow  field,  with  a  universal  model*  furthermore,  if  we 
try  to  construct  a  universal  model  by  using  one  flow  to  determine  one 
model  constant  and  another  flow  for  another  model  constant,  we  mav ,  to 
some  degree,  "contaminate"  the  model  and  render  it  invalid  for  either 
flow.  As  a  result,  the  model  mav  not  be  accurate  for  either  flow.  Ue 
will  show  below  that  this  "contamination"  is  stronger  than  has  been 
generally  recognized  for  two-equation  models. 

Since  a  complex  flow  field  contains  several  flow  zones  with  differ¬ 
ent  physical  structures,  and  since  universal  closure  models  are  usually 
"contaminated",  it  seems  logical  to  construct  models  In  which  each  flow 
zone  is  modeled  independently.  Here  we  define  a  flow  zone  as  a  distinct 
part  of  the  flow  field  with  definable  and  distinctive  physical  charac¬ 
teristics.  This  point  deserves  emphasis  —  the  models  are  tied  to  the 
local  cha  rac  t  e  r  i  st  i  cs  of  the  flow,  not  to  the  flow  as  a  whole  and  cer¬ 
tainly  not  to  a  wide  variety  of  flows.  Most  flows  contain  more  t han  one 
zone.  For  example,  a  diffuser  flow  with  separation  (Fig.  1.1)  can  be 
regarded  as  a  composite  of  six  flow  zones:  (A)  the  potential  core;  (B) 
the  attached  boundary  layer;  (C)  the  detachment  zone;  (D)  the  recir¬ 
culation  zone;  (F)  the  free  shear  layer;  (F)  the  reattachment  zone. 
With  zonal  models,  we  can  provide  a  separate  model  for  each  zone.  It 
should  be  much  easier  to  construct  accurate  models  for  zones  than  foi 
the  flow  field  as  a  whole,  and  the  models  can  be  much  simpler  than  a 
universal  model  would  need  to  be.  Equally  importantly,  there  is  little, 
if  any,  danger  of  contamination.  Fortunately,  there  is  a  limited  number 
of  types  of  flow  zones  in  flows  with  engineering  application.  There¬ 
fore,  the  task  of  zonal  modeling  is  large  but  finite.  We  believe  that 
the  zonal  modeling  approach  will  probably  yield  accurate  engineering 
methods  more  quickly  than  will  a  search  for  a  single  universal  closure 
model . 

In  this  approach,  the  definition  of  a  "zone"  must  he  subject  to 
pragmatic  testing.  At  the  beginning,  we  select  a  region  we  believe  has 
distinct  physics  as  a  candidate  zone.  We  then  construct  a  model  for 
that  zone.  If  the  model  works  within  the  desired  accuracy,  the  zone 
definition  is  accepted.  However,  the  process  mav  nor  end  here.  The 
following  possibilities  arise: 
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(a)  The  region  la  not  adequately  Modeled.  In  this  case,  the  zone 
May  nead  to  be  broken  Into  two  or  More  regions  each  with  a 
different  Model;  the  zonal  definition  is  then  narrowed. 

(b)  The  Model  works  not  only  for  the  region  intended,  but  for 
other  regions  as  well.  In  this  case,  the  definition  of  the 
zone  can  be  broadened. 

Thus  In  each  case,  we  are  guided  by  the  success  of  the  Model  and  not  by 

our  Initial  guesses.  In  the  end,  what  is  seen  as  a  zone  therefore  de¬ 

pends  strongly  on  the  Model  used. 

For  free  shear  flows,  we  began  by  adopting  the  classic  classifica¬ 
tion  acheMS,  that  Is,  jets,  wakes,  and  Mixing  layers  were  assumed  to  be 
distinct  zones.  Since  Jets  and  wakes  have  both  planar  and  axlsyMMetrlc 
realizations,  there  are  five  basic  cases.  In  each  case  there  is  a  near 
field  and  a  far  field.  This  suggests  that  as  many  as  ten  zones  night  be 
needed.  In  addition,  there  are  co-flowing  jets  and  nixing  layers  with 
variable  ratio  of  the  velocities  of  the  two  main  streans.  Clven  this 

range  of  cases.  It  was  not  entirely  clear  what  to  select  as  zones.  We 

started  by  assualng  the  need  for  a  separate  model  for  each  flow.  In  the 
end,  we  found  that  Many  cases  could  be  fitted  with  a  single  model,  so 
that  there  are  only  three  zones  in  the  final  model;  the  classification 
was  not  obvious  a  priori. 

This  process  sounds  awkward;  in  practice,  however,  it  not  only  led 
to  resolution  in  a  reasonable  tine,  but  also  was  instructive  in  provid¬ 
ing  Information  about  the  basic  parameters  needed  to  model  free  shear 
flows  and  hence  about  the  underlying  physics.  In  the  end,  we  were  able 
to  unify  the  Models  to  form  a  single  oodel  with  two  parameters.  This 
model  Is  more  successful  than  any  we  are  aware  of.  We  were  able  to  use 
linear  first-order  equations  for  all  the  readjustment  processes  (to  be 
described  later). 

If  the  flow  field  is  decomposed  into  zones,  ae  suggested  above, 
there  exist  regions  between  the  zones  which  can  be  called  "readjustment 
regions",  where  the  structure  adjusts  from  that  of  the  old  zone  to  that 
of  the  new  one.  Physically,  these  readjustments  are  Transitions  from 
one  type  of  structure  to  another.  The  modeling  of  these  readjustment 
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regions  must  reflect  this  changing  physics.  Pragmatism  requires  that, 
at  least  at  the  outset,  we  adopt  the  simplest  possible  model  for  the 

readjustment  regions.  It  also  means  that,  to  minimize  the  difficulty  of 
blending,  the  models  tor  various  zones  should  have  the  same  form;  for 

example,  they  may  all  be  of  the  K-e  type.  A  linear  blending  model 

using  a  first-order  lag  equation  meets  the  simplicity  criterion;  our 

expectation  Is  that  It  may  be  able  to  represent  most  readjustment 
regions  with  reasonable  accuracy.  Thus  far,  our  expectation  has  been 
fulfilled  as  will  be  reported  In  later  chapters.  However,  we  do  not 
rule  out  the  possibility  that  more  sophisticated  blending  models  may  be 
needed  in  some  cases. 

To  recapitulate,  in  constructing  zonal  models,  we  shall  be  guided 
bv  following  criteria: 

1.  ACCURACY .  Pach  model  should  represent  the  zone  It  Is  designed 
for  with  sufficient  accuracy.  Acceptable  accuracy  may  vary  with  appli¬ 
cation.  For  the  cases  of  Interest  to  the  authors,  5Z  accuracy  In  most 
Integral  parameters  is  a  reasonable  goal.  In  other  applications,  for 
example,  environmental  flows,  much  less  accuracy  may  suffice. 

-•  SIMPLICITY.  We  want  the  model  for  each  flow  zone  to  be  as  sim¬ 
ple  as  possible,  consistent  with  engineering  accuracy. 

1.  UNIFORMITY.  We  want  the  zonal  models  to  be  uniform  in  the 
sense  that  the  variables  used  In  each  zone  should  be  as  similar  as 
possible  in  order  to  simplify  the  zone  readjustments.  This  will  also 
help  to  make  model  development  and  assessment  relatively  straightfor¬ 
ward.  Two-equation  models  and,  par 1 1 cu lar 1 v ,  K-t  models,  will  he 
selected  as  the  basis  on  which  to  build  the  zonal  models;  this  will 
allow  maximum  advantage  to  be  taken  of  existing  models. 

a.  INDKPKNDKMCK.  We  want  the  ability  to  alter  the  model  In  one 
zone  without  affecting  the  models  in  other  zones.  This  will  allow 
r  inlng  for  each  zone  separately.  It  also  avoids  the  "contami  nafi  in" 
problem  discussed  above. 

‘.  n.OLt;  CMVNHTION  rn  PHYSICS.  As  mentioned  earlier,  turbulence 
models  should  reflect  the  physics  of  the  phenomena  t  hev  are  intended  t  > 
represent.  However,  incorporating  phvsics  into  a  model  is  difficult  and 


has  rarely  beer.  done.  Hiere  la  conalderable  accumulated  knowledge  on 
flow  structures  which  ahould  be  uaeful  In  thl8  effort. 

6.  NEGLIGIBLE  NUMERICAL  ERRORS.  This  Is  not  really  a  modeling 
Issue,  but  we  need  to  eliminate  uncertainty  due  to  the  numerical  errors 
In  coaputatlon.  This  la  necessary  If  one  Is  to  evaluate  the  models  ob¬ 
jectively  and  requires  the  construction  of  accurate  and  stable  nuaertcal 
aethods.  At  present,  adequate  nuaerlcal  aethods  exist  only  for  soae 
types  of  flows  and  s»re  should  become  available  In  the  next  few  years. 

1.2  An  Overview  of  Modeling  for  Turbulent  Free  Shear  Flows 


1.2.1  Background  for  Free  Shear  Flows 

Free  shear  flows  are  those  flows  in  which  there  is  no  direct  effect 
of  solid  boundaries.  This  class  of  flows  consists  of  five  basic  flows; 
the  aixing  layer,  and  the  plane  and  axlsymnetrlc  Jets  and  wakes.  As 
with  boundary  layers,  these  flows  have  a  single  predominant  flow  direc¬ 
tion  and  hence  the  shear  stresses  and  diffusion  fluxes  are  significant 
only  In  directions  perpendicular  to  the  predominant  direction.  Further¬ 
more,  pressure  gradients  normal  to  the  predominant  flow  direction  are 
usually  unimportant.  Consequently,  the  equations  governing  these  flows 
are  identical  to  those  of  boundary  layers.  For  this  reason,  the  free 
shear  flows  and  boundary  layers  are  sometimes  called  thin  shear  layers. 
The  free  shear  flows  are  all  Inherently  unstable  and  transition  from 
laminar  to  turbulent  occurs  at  Reynolds  numbers  between  15  and  30  based 
on  shear-layer  width.  It  Is  almost  impossible  to  maintain  a  laminar 
free  shear  flow  in  the  laboratory  and  all  free  shear  flows  of  practical 
engineering  interest  are  turbulent. 

Free  shear  flows  play  significant  roles  in  many  engineering  appli¬ 
cations.  Jet  engines.  Jet  ejectors,  fuel-oxidizer  mixing  in  propulsive 
devices,  and  wakes  behind  airplanes  and  submarines  are  a  few  examples. 
Also,  in  flows  with  separation,  free  shear  flow  may  represent  a  zone  of 
the  whole  flow.  Indeed,  free  shear  flow  zones  are  commonly  seen  in 
engineering  flows  (see  Kline  et  al.  (1981)). 

Turbulent  free  shear  flows  are  a  good  starting  point  for  studying 
turbulence  because  they  are  relatively  simple  and  hence  offer  a  better 
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chance  for  understanding.  Most  of  the  turbulence  structures  in  engin¬ 
eering  flows  are  not  well  understood  at  the  present  time.  When  extra 
strains,  such  as  the  effects  of  solid  walls  and  pressure  or  temperature 
gradients,  are  added  to  the  flows,  the  turbulence  evolution  becomes  more 
complicated  and  more  difficult  to  understand. 

In  classic  theory,  power  laws  describe  the  downstream  variation  of 
the  layer  width  and  the  velocity  decay  with  downstream  distance  for  many 
free  shear  flows  (see,  for  example,  Schlichting  (1979)  and  Tennek.es  and 
Lumley  (1972)).  These  power  laws  are  exact  for  laminar  flows  and  are 
based  on  dimensional  arguments  and  the  momentum  integral  equation  for 
turbulent  flows.  Although  not  exact,  power  laws  do  provide  a  guideline 
of  the  behavior  of  the  flows  and  experimentalists  often  use  power  laws 
for  fitting  data. 

With  present  computers,  we  cannot  solve  the  full  time-dependent 
Navler-Stokes  equations  (full  turbulence  simulation),  which  are  the  only 
truly  unlvenal  "model"  for  all  the  turbulent  flows,  except  for  very 
simple  flows.  Fortunately,  the  details  of  a  turbulent  flow  are  rarely 
needed;  time-averaged  quantities  usually  Buffice  even  when  the  mean  flow 
Is  unsteady.  Time-averaged  equations  are  therefore  used  for  most  engin¬ 
eering  calculations.  In  the  process  of  time  averaging,  however,  most  of 
the  detailed  Information  contained  in  the  Navier-Stokes  equations  Is 
lost  and,  as  a  result  of  the  nonlinearity,  the  time-averaged  equations 
do  not  form  a  closed  set  of  equations.  Closure  of  the  set  of  time- 
averaged  equations  requires  a  turbulence  model  describing  the  Reynolds 
stresses  In  terms  of  quantities  that  can  be  calculated. 


1.2.2  Eddy  Viscosity  Model 

Boussinesq  (1877)  Initiated  the  concept  of  eddy  viscosity;  It  is 
still  widely  used  In  turbulence  modeling.  He  suggested  that  the  effec¬ 
tive  turbulent  shear  stress,  arising  from  the  cross-correlation  of  the 
fluctuating  velocities,  be  replaced  by  the  product  of  mean  velocity 
gradient  and  a  quantity  called  the  "eddy  viscosity"  or  "turbulent  vis¬ 
cosity".  This  model  resembles  the  stress-strain  relationship  In  laminar 
flows.  Unlike  the  molecular  viscosity,  the  eddy  viscosity  Is  not  a 
property  of  the  fluid.  Its  value  varies  from  point  to  point  In  the 


flow,  and  la  largely  deterained  by  Che  characteristics  of  the  turbulence 
at  the  point  In  question.  The  Introduction  of  the  eddy  viscosity 
provides  a  framework  for  turbulence  modeling,  but  It  does  not  itself 
constitute  a  model;  there  remains  the  task  of  expressing  the  eddy 
viscosity  in  terms  of  known  or  calculable  quantities. 

The  mixing  length  hypothesis,  introduced  by  Prandtl  (1925),  has 
been  a  major  success  in  turbulence  modeling.  It  assumes  that  the  eddy 
viscosity  is  proportional  to  the  product  of  the  density,  a  length  scale, 
and  a  velocity  scale.  For  two-dimensional  flows,  it  further  assumes 
that  the  velocity  scale  can  be  replaced  by  the  product  of  the  length 
scale  and  the  mean  velocity  gradient.  This  reduces  the  number  of  un¬ 
knowns  to  one,  a  length  scale  known  as  the  mixing  length.  For  flows 
bounded  by  walls,  the  length  scale  Is  proportional  to  the  distance  from 
the  wall  in  regions  very  close  to  the  wall  (inner  layer)  and  to  the 
layer  thickness  where  the  wall  effect  is  small  (outer  layer).  For  tur¬ 
bulent  free  shear  flows,  the  length  scale  should  be  proportional  to  the 
width  of  the  layer,  but  the  proportionality  constant  varies  from  flow  to 
flow  (see,  for  example.  Launder  and  Spalding  (1972)).  There  is  no  way 
of  predicting  the  constant  for  a  flow  that  has  not  been  measured.  This 
lack  of  generality  is  a  major  drawback  of  this  method. 

For  free  shear  flows,  the  mixing  length  model  assumes  constant 
length  scale  across  the  layer  at  a  given  downstream  location.  Another 
class  of  models  assumes  the  eddy  viscosity  constant  across  the  layer. 
In  these  models,  the  velocity  scale  is  taken  as  the  maximum  velocity 
difference  in  the  layer  rather  than  the  product  of  length  scale  and 
velocity  gradient.  Generally,  the  constant  eddy  viscosity  model  per¬ 
forms  less  satisfactorily  than  the  mixing  length  model.  However,  the 
constant  eddy  viscosity  model  allows  similarity  solutions  for  which  the 
problem  can  be  reduced  to  an  ordinary  differential  equation  (Schlichting 
(1979)  and  White  (1974)).  This  was  a  major  advantage  when  high- 
performance  computers  and  accurate  numerical  methods  for  solving  partial 
differential  equations  were  not  available. 


1.2.3  Two-Equation  Models 


It  was  not  until  the  late  1960s  that  differential  equations  were 
used  to  calculate  the  velocity  and  length  scales  that  determine  the  eddy 
viscosity.  In  these  models,  the  velocity  scale  is  proportional  to  the 
square  root  of  the  turbulent  kinetic  energy  which  is  obtained  by  solving 
a  differential  equation.  If  the  length  scale  is  prescribed,  we  have  the 
so-called  one-equation  models.  Alternatively,  the  length  scale  may  be 
derived  from  another  differential  equation;  we  then  have  two-equation 
models.  Several  length-scale-determining  equations  have  been  proposed. 
The  one  which  describes  the  turbulent  energy  dissipation  rate  has  been 
favored  by  many  researchers.  The  reason  lies  in  the  fact  that  this 
quantity  appears  in  the  turbulent  kinetic  energy  equation.  This  model 
is  called  the  K-e  model  and  has  enjoyed  a  great  deal  of  success  since 
Jones  and  Launder  (1972)  first  proposed  it.  Rodi  (1972)  has  used  it 
extensively  in  the  prediction  of  free  shear  flows.  Two-equation  models 
do  not  require  flow-specific  length  scale  information  and  offer  the  pos¬ 
sibility  of  achieving  the  generality  that  the  mixing  length  and  constant 
eddy  viscosity  models  lack.  However,  the  performance  of  the  K-e  model 
falls  short  of  the  desired  accuracy.  For  example,  the  standard  K-e 
model  over-predicts  the  spreading  rate  of  an  axisymmetric  jet  by  30-35% 
and  under-predicts  the  plane  wake  spreading  rate  by  nearly  40%.  The 
reason  for  the  poor  predictions  of  these  flows  may  be  contamination.  On 
the  other  hand,  only  slight  adjustments  are  required  to  improve  the 
model.  Details  will  be  given  in  later  chapters. 

This  provides  the  motivation  for  the  present  study  of  free  shear 
flows.  We  will  apply  the  zonal  modeling  strategy  to  the  development  of 
improved  models.  We  begin  with  a  close  examination  of  the  physics  of 
individual  zones  and  then  incorporate  as  much  of  what  we  learn  as  pos¬ 
sible  into  the  model.  The  standard  K-e  model  is  adopted  as  a  base 
model  since  it  predicts  many  flows  reasonably  well.  Systematic  modifi¬ 
cation  of  the  model  constants  for  well-defined  zones  is  performed  and 
then  the  zonal  models  are  blended  to  allow  readjustments  between  the 
zones.  It  is  found  that  one  primary  and  one  secondary  parameter  suffice 
to  classify  and  predict  all  free  shear  flows  considered.  The  primary 
parameter  is  the  free  stream  velocity  divided  by  the  maximum  velocity 
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difference  in  the  layer;  it  represents  the  ratio  of  the  inertial  force 
to  the  driving  force.  This  parameter  is  zero  for  "pure"  jets.  It 
varies  from  a  small  value  to  Infinity  in  co-flowing  jets  or  wakes.  In 
mixing  layers,  this  parameter  is  fixed  at  some  finite  value,  determined 
by  the  velocity  ratio  of  the  two  streams,  and  can  take  on  values  from 
one  to  infinity.  As  this  parameter  increases,  the  driving  force  becomes 
smaller  compared  to  the  inertial  force,  it  is  harder  for  eddies  to  pene¬ 
trate  into  the  free  stream,  and  the  flow  becomes  more  parallel.  In  this 
limit,  the  dominant  mechanism  of  the  turbulence  transport  is  diffusion. 
This  effect  is  simulated  by  changing  two  model  constants. 

The  secondary  parameter  is  a  measure  of  lateral  vortex  stretching. 
This  vortex  stretching  reduces  the  correlation  between  the  shear  stress 
and  the  turbulent  kinetic  energy.  For  example,  the  maximum  value  of 
lu'v'l/K  is  only  about  0.23  in  axisymmetric  "pure"  jet,  in  which  vortex 
stretching  is  strongest,  while  it  is  about  0.3  in  all  plane  flows.  This 
effect  disappears  in  far-downstream  axisymmetric  flows.  Change  of  one 
model  constant  is  required  to  incorporate  this  effect. 

Experiments  are  the  basis  of  our  understanding  of  turbulence  phe¬ 
nomena  and  provide  the  data  required  to  fix  the  model  constants  or  func¬ 
tions.  Therefore,  they  are  crucial  in  the  development  of  zonal  models 
and  readjustments.  There  are  extensive  experiments  on  turbulent  free 
shear  flows.  They  will  be  carefully  reviewed  below. 

In  Chapter  II,  the  standard  "universal"  model  is  reviewed  and  the 
zonal  modeling  concept  is  tested  on  three  homogeneous  flows.  In  Chapter 
III,  a  survey  of  free  shear  flow  experiments  is  performed.  In  Chapters 
IV  and  V,  a  zonal  modeling  approach  is  used  for  the  development  of  the 
model  for  plane  and  axisymmetric  free  shear  flows  respectively.  In 
Chapter  VI,  a  summary  of  the  new  model  that  brings  together  the  plane 
and  axisymmetric  flows  is  presented. 
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Chapter  II 

ZONAL  MODELING  OF  HOMOGENEOUS  FLOWS 


8 


Cv 

$ 


2.1  Background 

Since  the  introduction  of  the  idea  of  zonal  modeling  (Kline  et  al. 
(1981)),  many  researchers  have  shown  interest  in  the  approach.  Others 
have  criticized  it  and  raised  questions.  It  was  not  clear  whether  this 
concept  would  work,  although  it  appeared  to  have  promise.  It  was  there¬ 
fore  decided  to  make  a  quick  test  of  the  idea  before  initiating  a  sig¬ 
nificant  effort. 

Homogeneous  flows  were  selected  for  this  task.  Homogeneous  flows 
are  ones  in  which  statistical  turbulence  quantities  are  independent  of 
spatial  position  and,  therefore,  functions  of  time  only;  they  may  con¬ 
tain  a  mean  flow  with  a  constant  rate  of  strain.  Since  nearly  all 
engineering  flows  are  inhomogeneous,  homogeneous  flows  have  little 
direct  engineering  application.  Nevertheless,  they  are  simple  turbulent 
flows  and,  consequently,  it  is  easier  to  construct  models  for  them  than 
for  inhomogeneous  ones.  We  shall  use  these  homogeneous  flows,  because 
they  provide  an  easy  test  of  the  zonal  modeling  concept  and  will  lay  a 
foundation  for  later  development  of  zonal  modeling  for  more  complex 
flows . 

In  this  chapter,  zonal  models  for  three  homogeneous  turbulent  flows 
are  developed.  These  are:  homogeneous  shear;  plane  strain;  and  axisym- 
metric  strain  flows.  All  of  these  flows  start  with  Isotropic  turbulence 
as  the  initial  condition.  When  shear  or  strain  is  applied  to  the  flow, 
the  turbulence  begins  to  depart  from  isotropy  and  undergoes  a  transition 
or  readjustment  to  a  state  which  depends  on  the  nature  of  the  applied 
strain.  After  a  long  enough  time,  the  turbulence  reaches  an  asymptotic 
state.  These  flows  can  be  regarded  as  consisting  of  two  zones  separated 
by  a  readjustment  stage;  the  two  zones  are  the  initial  isotropic  flow 
zone  and  the  fully-developed  zone.  Since  homogeneous  flows  are  time- 
developing  flows,  these  flow  zones  occur  in  time  rather  than  in  space. 
Although  the  dimension  of  variation  is  different  from  the  flow  zones 
described  in  Chapter  I,  the  concept  is  much  the  same. 
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2.2  The  Standard  K-e  Model 


Before  developing  the  new  zonal  models  for  homogeneous  flows,  it  is 
worthwhile  to  examine  the  existing  "universal"  K-e  model  (see,  for 
example,  Rodi,  1980).  In  this  model,  the  Reynolds  stress  tensor  is  mod¬ 
eled  using  the  eddy  viscosity  concept: 


3U  3U 


-u  u 

1  j 


/  i  j  \  2 

- +  — ~  ~  r  6  K 

t  \  3x.  3 x  /  3  ij 


(2.2.1) 


The  eddy  viscosity  is  given  by: 


(2.2.2) 


The  turbulent  kinetic  energy,  K,  and  rate  of  dissipation  of  the  turbu¬ 
lent  kinetic  energy,  e,  are  governed  by  the  differential  equations: 


(2.2.3) 


where  Uj  is  the  mean  velocity  in  the  i-direction.  C  ,  C  ,  C  ,  a,  , 

U  £  1  K 

and  ae  are  model  constants.  The  first  two  terms  on  the  rig ltt  hand 
side  of  each  model  equation  represent  production  and  destruction  of  the 
quantity  whose  evolution  the  equation  describes.  The  last  terms  of 
these  two  equations  represent  diffusion,  which  plays  no  role  in  homo¬ 
geneous  f lows . 

Values  of  C  *  0.09,  C  =  1.44,  C_  =  1.92,  a*.  “  1.0,  and 

"  E2 

and  oe  “  1.3  are  widely  used  for  the  model  constants.  For  reference, 
we  shall  call  these  the  "standard  values"  and  the  model  equations  with 
these  values  the  "standard  K-e  model". 


The  five  standard  values  of  model  constants  were  determined  as  fol 

lows: 


are  zero  so 


In  homogeneous  grid  turbulence,  diffusion  and  production 

C_  is  the  only  constant  that  plays  a  role  in  Eqs.  (2.2.3) 
e2 


and  (2.2.4).  From  the  measured  rate  of  decay  of  turbulent  kinetic 
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energy,  C£^  can  be  determined  and  is  found  to  lie  in  the  range  1.8  to 
2.0;  as  noted  above,  1.92  is  the  most  commonly  used  value. 

2.  Cy:  In  local-equilibrium  shear  layers,  production  and  dis¬ 

sipation  are  essentially  equal,  and  Eqs.  (2.2.1)  and  (2.2.2)  can  be  com- 

-  2  - 

bined  to  yield  C, ,  =  fu  u  /K)  .  Measurements  show  that  -u  u  /K  «0.3 
l1  M2'  12 

in  these  flows,  so  that  Cy  =  0.09. 

3.  C£  ^ ;  In  near-wall  regions  of  boundary  layers,  the  velocity 
profile  is  nearly  logarithmic,  production  is  approximately  equal  to 
dissipation,  and  the  convection  of  dissipation  is  negligible.  In  this 
situation,  Eq.  (2.2.4)  reduces  to 


a  /C 

e  u 


where  k  =  0.41. 


4.  Finally,  the  diffusion  constants  and  a£  were  assumed  to 
be  close  to  unity  and  they  were  tuned  by  computer  optimization.  C£^ 
is  then  obtained  from  the  eequation  above. 

Note  that  at  least  three  different  types  of  flows  (homogeneous  grid 
turbulence,  equilibrium  free  shear  layer,  and  boundary  layer)  were  used 
to  determine  the  model  constants.  This  use  of  different  flows  causes 
the  type  of  contamination  described  in  the  Introduction.  The  standard 
K-e  model  probably  predicts  boundary-layer  flows  better  than  other 
flows  because  two  of  the  model  constants  were  set  using  data  from 
boundary-layer  flows.  However,  the  standard  model,  although  widely 
used,  cannot  be  expected  to  work  well  "universally".  We  shall  assess 
how  well  (or  badly)  it  does  for  a  number  of  flows  and  several  variables 
in  this  work. 

Applied  to  homogeneous  flows,  Eqs.  (2.2.3)  and  (2.2.4)  reduce  to 
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and  only  three  model  constants,  C^,  Ce^,  and  C^.  reBain  *n  the 

equations.  Note  that  the  physics  from  which  C  and  C*  were  deter- 

V 

mined  are  not  appropriate  to  homogeneous  flows.  Therefore  it  is  ques¬ 
tionable  whether  the  standard  K-e  model  is  acccurate  for  such  flows. 
In  fact,  its  performance  for  the  homogeneous  flow  cases  in  the  1980-81 
Stanford  Conference  was  not  good.  Some  examples  of  the  performance  of 
this  model  for  homogeneous  flows  will  be  given  later. 

In  the  present  study,  zonal  models  for  homogeneous  shear  flow, 
homogeneous  axisymmetrlc  strain  flow,  and  homogeneous  plane  strain  flow 
are  constructed  as  described  below.  The  inputs  used  to  build  the  models 
came  from  the  full  turbulence  simulations  by  Rogallo  (1981)  and  Lee  and 
Reynolds  (1985).  These  zonal  models  perform  considerably  better  than 
the  standard  K-e  model  in  predictions  of  these  flows,  since  they  have 
been  tuned  for  these  xows. 


2.3  Zonal  Models  for  Homogeneous  Flows 


2.3.1  Homogeneous  Shear  Flow 


Homogeneous  sheared  turbulence  is  a  flow  in  which  Initially  iso¬ 
tropic  turbulence  is  subjected  to  a  mean  shear  flow.  After  application 
of  the  shear,  the  structure  of  the  turbulence  becomes  anisotropic.  The 
flow  gradually  adjusts  to  the  shear  and,  after  a  long  enough  time, 
appears  to  reach  a  fully-developed  or  asymptotic  state  in  which  the 
physics  no  longer  changes.  This  flow  can  be  regarded  as  consisting  of 
two  flow  zones  connected  by  a  readjustment  stage.  The  two  zones  are  the 
isotropic  initial  flow  zone  and  the  fully-developed  shear  flow  zone. 
These  two  flow  zones  occur  in  the  time  domain.  Between  these  two  time 
periods  is  the  readjustment  stage  where  the  change  of  physics  takes 
place.  The  readjustment  Is  the  most  interesting  part  of  this  flow  and 
the  part  to  which  most  of  our  attention  will  be  devoted. 


14 


where 

stress 


r  -  dU/dy 

s 

component  is 


is  the  shear  rate.  The  most  important  Reynolds 
U1U2*  Initially,  in  the  isotropic  state. 


v3 


2 

3  6ijK 


(2.3.2) 


i.e.,  uiu2  “  0»  and  in  the  fully-developed  state, 

u  u  **  -0.32K  (2.3.3) 

1  2 

T  ~2  ~~ 2. 

where  K  =  (u^  +  +  u^J/2  is  the  turbulent  kinetic  energy.  Equation 

(2.3.3)  correlates  the  experimental  data  extremely  well.  No  modeling  is 
needed  for  the  other  components  of  the  Reynolds  stress.  We  shall  there¬ 
fore  adopt  this  as  our  model  for  the  fully-developed  state. 

In  the  readjustment  region,  we  shall  assume  that  u^u^  can 
represented  by  a  simple  blend  of  the  initial  and  final  states: 

=  -0.32  X  K  (2.3.4) 

where  X  is  the  readjustment  parameter  utilized  to  model  the  change  of 
physics.  It  is  defined  to  have  the  range  0  _>.  X  _>  1. 

In  line  with  the  philosophy  of  simplicity  expounded  in  the  intro¬ 
duction,  we  shall  assume  the  readjustment  parameter  obeys  a  first-order 
lag  equation,  specifically,  X  is  determined  by  the  first-order  ordi¬ 
nary  differential  equation: 


—  -  b  T  (1  -  X)  (2.3.5) 

dt  s  s 

where  bg  is  a  dimensionless  constant.  This  choice  of  lag  equation 
makes  the  rate  of  change  of  the  physics  greatest  for  small  X.  This 
behavior  accords  with  observation.  The  solution  to  Eq.  (2.3.5)  is 

-b  r  t 

X  «  1  -  e  3  8  (2.3.6) 


By  fitting  the  experimental  data  from  Case  376A  of  the  1980-81  confer¬ 
ence,  it  is  found  that 


The  model  equation  for  K  can  be  written,  with  the  help  of  Eq.  (2.3.4), 
as 

-  0.32  XT  K  -  e  (2.3.8) 

at  s 

The  first  term  on  the  right  hand  side  is  the  rate  of  production  of  kin¬ 
etic  energy.  When  X  -  0,  the  production  term  vanishes  and  Eq.  (2.3.8) 
reduces  to  a  form  appropriate  to  isotropic  flows.  When  X  -  1,  Eq . 
(2.3.8)  models  fully-developed  shear  flows. 

Following  the  form  used  in  the  standard  K-e  model,  the  e  model 
equation  can  be  written  as 

2 

—  -  C  0.32  xr  e  -  C  —  (2.3.9) 

dt  e  s  e  K 

1  2 

Fitting  Rogallo's  (1981)  full  simulation,  we  find 


1.406 


(2.3. 10) 


and 


n_ 

6 


(2.3.  11  ) 


This  value  of  has  some  analytical  underpinning  and  Is  favored  by 
Reynolds  (1976).  When  X  -  0,  Eq.  (2.3.9)  reduces  to  the  model  for  the 
decay  of  Isotropic  turbulence  (see  Reynolds  (1976  )),  and  as  X  Is  In¬ 
creased  from  zero,  the  production  of  dissipation  term  (the  first  term) 
Increases . 

The  Reynolds  stress  model,  Eq.  (2.3.4)  and  the  K  and  e  equa¬ 
tions  (2.3.8)  and  (2.3.9),  together  with  the  equation  for  the  readjust¬ 
ment  parameter,  Eq.  (2.3.6)  constitute  the  zonal  model  for  homogeneous 
shear  flows  with  isotropic  initial  fields.  When  X  -  0,  the  system 
reduces  to  the  model  for  isotropic  flow  and,  when  X  ■  I,  to  the  fully- 
developed  shear  flow  model.  As  X  moves  between  zero  and  unity,  the 
system  undergoes  a  transition  from  one  form  to  the  other  and  allows  an 
adequate  readjustment  from  one  zone  to  the  other.  Results  obtained  with 
this  model  will  be  given  below. 


2.3.2  Homogeneous  Axl symmetric  Strain  Flow 


Again,  In  this  flow  the  Initial  velocity  field  Is  Isotropic  turbu¬ 
lence.  At  time  t  “  it  begins  to  be  stretched  in  one  direction  and 

compressed  equally  in  the  other  two  directions.  This  axl svmaet r  i  c 
strain  causes  anisotropy;  turbulence  fluctuations  in  the  stretching 
direction  are  reduced  while  those  in  the  compression  directions  are 
Increased.  This  is  a  result  of  the  vortex  filaments  aligning  with  the 
stretching  direction.  After  the  strain  is  applied  for  a  sufficient  lv 
long  rime,  the  Reynolds  stress  reaches  a  fully-developed  structure. 
Between  the  initial  and  asymptotic  states,  the  Reynolds  stress  undergoes 
a  transition  from  one  tvpe  of  structure  to  the  other.  Again  the  time 
historv  can  be  decomposed  into  two  zones  separated  by  a  readjustment 
s  t  age . 

The  strain  tensor  for  axi symmet r i c  strain  is 


where  r  -  dt'  dx  is  the  strain  rate.  The  Revnolds  stress  in  the  iso¬ 
tropic  state  is  given  bv  F.q .  '.'.1.2).  After  the  turbulence  undergoes 
axl symmet  r i c  strain  for  i  long  period,  the  Revnolds  stress  tensor 
reaches  an  asvmptotic  state  In  which 


a 


which  shows  that  the  component  energy  approaches  zero  in  the  stretching 
direction  an'  l'  :n  ear*-  >t  'he  compression  directions.  Again  assiming 
a  monot  >n  i ,  ■  vui  r  i  >n  *  t  m  the  initial  state  des  r  i  bed  bv  K.q . 
and  ’he  state  ie'ined  bv  Kq .  3  .  I .  H  ' ,  we  have 
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The  behavior  >f  *,  the  readjustment  parameter,  Is  again  assumed 
to  be  described  bv  a  first -order  lag  equation 
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where  b  (  Is  a  constant  which  need  nor  be  the  saise  as  of  Kq . 

(  2 .  1 .  S  ) .  The  solution  to  Rq  .  (  2  .  1 .  1  *>  )  Is 
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Fitting  Roga  1  1  o  '  s  9  8  j  '  til.  stmu 1  at  :  in,  we  ‘  :  nd  that 


With  R.q  .  1. Is  > ,  t  he  •’  mode  1  »»,,  c  l  on  t  s 
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Again  we  ise  t  he  standard  ’  >  r  m  of  the  >  model  equation  for  this  flow 
and  a  r  r  1  v a  t 
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However,  within  the  philosophy  if  /onal  modeling,  there  is  no  reason  whv 
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2.1.3  Homogeneous  Plane  Strain  Flow 


In  this  flow,  the  initially  isotropic  flow  field  is  stretched  in 
one  direction  and  compressed  in  a  second  direction,  with  the  third 
direction  being  neither  stretched  nor  compressed.  The  fluctuations 
diminish  In  the  stretching  direction  and  increase  in  the  compression 
direction  while  remaining  close  to  2  K / 3  in  the  neutral  direction. 
F.ventuallv,  the  Reynolds  stress  approaches  a  fully-developed  structure. 
Thus,  th  two-rone  concept  In  the  time  domain  applies  to  this  flow. 

The  strain  tensor  for  plane  strain  is 
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where  “  dW 'dz  -  -dV 'dv  is  the  strain  rate.  We  .-an  model  the  n  >rrta 

Reynolds  stresses  for  the  component  energies')  bv 
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where  a  |  and  an  are  constants,  and  0  s  A  s  1  is  the  readjustment 
parameter.  ^.s  in  the  above  two  flows,  A  Is  found  from  a  first-order 
linear  ordinary  differential  equation: 
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with  the  solution 
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The  constant  b  may  be  different  than  in  the  other  flows.  Note  that, 
P 

In  thla  case,  unlike  the  previous  ones,  for  fixed  total  strain  r  t  , 

~~2  ~ 2  2 

Uj  and  Uj  vary  with  the  dimens lonless  strain  rate  t  while  u^ 

is  Independent  of  the  strain  rate. 

To  determine  the  three  constants,  aj,  a2>  and  b  In  the  isodel, 
we  used  the  following  procedure: 

2 

1.  By  fitting  the  data  of  iw'K  fro*  l.ee  and  Reynolds'  (1985) 
full  simulation.  It  is  found  that 


b  •  l  .  1 
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(2.3. 28  ) 


2.  From  rapid  distortion  theory  for  plane  strain,  Batchelor  and 

~ 2  ~2 

Proudaan  (1954),  both  u^  and  u^  tend  to  approach  K  for  large  values 
of  r  t  when  Is  large.  This  is  because,  under  these  conditions, 

the  vortex  filaments  align  with  the  stretching  direction  so  rapidly  that 
the  component  energies  In  the  other  two  directions  each  share  half  of 
the  total  energy  2K.  This  yields: 

a,  -  3  v2.3.29) 

2 

3.  Finally,  by  fitting  Lee  and  Reynolds'  (1985)  data  for  u ? /K , 
the  component  energy  In  the  compression  direction,  it  is  found  that 


a  ,  *  26 .  ) 


(  2  .  3 .  30  ) 


With  the  developed  normal  stress  functions,  F.qs.  (2.3.23),  (2.3.24),  and 
(2.3.25),  the  two  model  equations  for  K  and  t  can  be  written  as 
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(2.  3.  32  ) 


Fitting  the  full  simulation  data  of  Lee  and  Reynolds  (1985),  we  obtain 


C  -  1.75  (2.3.33) 

‘-l 

and 

C  -  —  (2.3.34) 

'’2  6 

Therefore,  the  zonal  model  for  the  plane  strain  consists  of  the  model 
equations  (2.1.31)  and  (2.3.32)  the  Reynolds  normal  stress  functions, 
F.qs.  12.3.23),  '2.1.24),  and  (2.3.25),  and  the  readjustment  parameter, 

Fq.  (2.1.27). 


^.1.4  Summa  r  v 


In  summarv,  the  modal  constants  b,  C  ,  and  C  for  the  three 

'  1  L  2 

zonal  models  developed  here  are  tabulated  as  below.  The  constants  C 
and  f  In  the  standard  K-r  model  are  also  Included  for  comparison. 
Note  that  the  > -model  equations  (2.2.6),  (2.3.9),  (2.3.19)  and  (2.3.32) 
ar>*  all  written  in  a  consistent  form: 
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d  K 
dt 
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where  P 
energy . 


-  V,  s i j 


is  the  production  rate  of  turbulent  kinetic 
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Comparison  of  Model  Constants  (Homogeneous  Flows) 


Standard 

Zonal  Model 

k-r  Model 

Homogeneous 

Shear 

FI  ow 

Homogeneous 
4xlsymraet  rlc 
Strain 

Homogeneous 
p  1  ane 

St  ra 1 n 

Lag  constant 
b 

1  .  5 

1.14 

1  .  1 

r 

‘  l 

1  .  44 

1  .  4 Of) 
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1 

1  .  7  5 
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Note  that,  although  our  models  are  simple,  they  apply  strictly  only  to 
cases  In  which  one  fixed  type  of  strain  is  applied  to  initially  iso¬ 
tropic  turbulence.  More  complex  cases,  such  as  those  with  initially 
anisotropic  states  or  strains  which  vary  with  time  may  not  be  amenable 
to  treatment  by  these  models.  However,  as  our  purpose  is  to  demonstrate 
the  inadequacy  of  the  standard  K-e  model  and  the  validity  of  the  zonal 
concept  and  not  to  develop  a  general  model  of  homogeneous  flows,  these 
are  not  serious  deficiencies  for  our  purposes.  Certainly,  a  universal 
model  would  be  valuable  but,  as  noted  earlier,  it  would  also  need  to  be 
much  more  complicated.  This  is  the  trade-off  we  intended  to  make;  we 
give  up  hope  of  universality  for  simplicity  while  maintaining  accuracy. 

2.4  Results  and  Comparisons 

Due  to  the  nature  of  homogeneous  flows,  there  is  no  convection  or 
diffusion  of  turbulence.  This  fact  allows  reduction  of  the  model  equa¬ 
tions  from  their  normal  partial  differential  equation  form  to  a  much 
simpler  ordinary  differential  equation  form.  These  model  equations  can 
be  easily  solved  by  numerical  Integration. 

The  Runge-Kutta  fourth-order  integration  scheme  seems  ideal  for 
this  integration,  and  is  used  In  the  present  study.  This  scheme  is  of 
' igh  accuracy  so  that  numerical  errors  are  eliminated  and  the  true  per¬ 
formance  of  the  model  can  be  evaluated.  The  three  zonal  models  for 
homogeneous  flows  described  above  are  tested.  The  results  obtained  are 
compared  with  those  predicted  by  the  standard  K-e  model.  Test  cases 
are  provided  by  Rogallo's  (1981)  and  Lee  and  Reynolds'  (1985)  full  simu¬ 
lations,  and  the  experiments  accepted  for  the  1Q80-H1  Stanford  Confer¬ 
ence  '"Kline  et  a  1 .  (1981)). 

2.4.1  Homogeneous  Shear  Flows 

1  .  R.  Ill's  Full  Simulation 

Duiia!  lo's  (198!)  full  simulation  of  homogeneous  shear  flows  were 
used  to  test  both  the  zonal  model  and  the  standard  K-e  model.  The 
re-u’ts  are  shown  In  Ffzs.  2. 1-2.4  for  various  shear  rates  and  viscosi¬ 
ties.  The  initial  turbulence  of  the  full  simulation  is  unrealistic 
(undev.-l  iped  isotropl'  fie’d,  square  spectrum),  and  it  takes  time  for 


w 


w 


the  simulation  to  develop  a  realization  of  real  turbulence.  Therefore, 
the  prediction  of  these  flows  by  both  models  was  begun  at  the  third 
time-step  of  Rogallo's  data.  Also,  the  last  data  points  of  the  full 
simulations  are  unreliable.  This  is  because  the  eddy  sizes  grow  with 
time,  and  when  they  approach  the  dimension  of  the  box  used  in  the 
simulation,  the  periodic  boundary  conditions  used  by  Rogallo  become 
invalid.  For  each  case,  twice  the  kinetic  energy  (i.e.,  2K)  and 
Reynolds  stress  -u^u^  vs*  c*me  are  plotted  in  Fig.  (a)  and  dissipation 
level  e  vs.  time  in  Fig.  (b).  For  all  cases,  the  zonal  model  predic¬ 


tions  agree  satisfactorily  with  the  full  simulation  data.  On  the  other 
hand,  the  standard  K-e  model  predictions  of  the  three  quantities  all 
grow  too  fast.  The  assumed  functional  form  of  the  Reynolds  stress  for 
the  K-e  model,  i.e.. 


-u  u 
1  2 


2 

K  dU 
U  e  dy 


is  responsible  for  the  overprediction  of  the  growth  by  the  K-e  model, 
because,  at  the  given  levels  of  K,  e,  and  dU/dy,  this  Reynolds  stress 
function  gives 


-u  u  >  0.32  K 

1  2 


which  makes  the  production  rate  too  high  and  hence  overpredicts  all  the 
three  quantities. 


2.  1980-81  Conference  Cases  376A  and  376B 


The  zonal  model  and  the  standard  K-e  model  were  both  employed  to 
predict  the  two  homogeneous  shear  flow  cases  (376A,  376B)  from  the  1980- 
81  Stanford  Conference.  The  results  are  shown  in  Figs.  2.5  and  2.6, 
respectively.  The  zonal  model  agrees  quite  well  with  the  experimental 
data  for  both  cases.  The  standard  K-e  model  underest imates  K  and 


-  for  case  376A  and  overestimates  them  for  case  376B.  The 


Reynolds  stress  function  again  sepms  to  be  the  source  of  the  under/over 


predictions  by  the  K-e  model;  the  predicted  values  of  -uiu2  are 


considerably  smaller  and  greater  than  0.32K  for  cases  376A  and  376B, 


respectively.  This  violates  the  established  result  that  -  u.u0 
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0.32K.  In  the  zonal  model  for  the  shear  flows,  however,  -u^u2  “  0*32K 
Is  explicitly  built  into  the  Reynolds  stress  function. 

2.4.2  Homogeneous  Axisymmetric  Strain  Flows 

Eight  cases  of  Rogallo's  (1981)  full  simulations  of  homogeneous 
axisymmetric  strain  flows  with  various  strain  rates  and  two  different 
viscosities  were  used  for  testing  the  zonal  model  and  the  standard 
K-e  model.  Figures  2.7-2.14  compare  the  predictions  of  both  models 
with  the  full  simulation.  For  each  case,  Fig.  (a)  shows  the  total  and 
component  energies  of  the  turbulence  vs.  time  and  (b)  the  dissipation 
rate.  The  zonal  model  performs  quite  well  for  the  whole  range  of  strain 
rates.  For  the  low-strain-rate  cases,  the  standard  K-e  model  accu¬ 
rately  predicts  K  but  not  the  component  energies.  The  dissipation 
level  is  also  incorrect  since  it  must  compensate  for  the  erroneous  pro¬ 
duction  rate.  When  the  strain  rate  is  high,  the  K-e  model  can  no 
longer  provide  accurate  prediction  of  K.  It  is  obvious  that  the 
assumed  Reynolds  stress  function  in  the  K-e  model  is  incorrect  for 
these  flows. 

No  appropriate  experiment  for  this  type  of  flows  is  available  from 
the  1980-81  Stanford  Conference  to  test  the  zonal  and  K-e  models. 

2.4.3  Homogeneous  Plane  Strain  Flows 

1 .  Lee  and  Reynolds'  Full  Simulation 

Lee  and  Reynolds'  (1985)  simulation  provides  six  cases  of  plane 
strain  flows  with  the  strain  rate  varying  from  0.65  sec  *  to  100  sec  *. 
The  zonal  and  standard  K-e  models  are  tested  against  these  data.  The 
results  of  total  and  component  energies  and  dissipation  vs.  time  are 
presented  in  Figs.  2.15-2.20.  The  zonal  model  predicts  the  flows  accu¬ 
rately  over  the  entire  range  of  strain  rates  simulated.  The  K-e 
model,  however,  does  not  perform  well  in  any  case,  especially  for  high 
strain  rates;  the  component  energies  are  all  seriously  in  error.  When 
strain  rate  is  high,  the  smallest  component  energy  (in  the  stretching 
direction)  even  becomes  negative,  which  is  physically  impossible.  This 
again  suggests  that  the  K-e  model  Reynolds  stress  function  is  not 
suitable  for  plane  straining  flows  and  further  reinforces  the  remark  on 
"contamination”  in  the  introduction. 


2.  1980-81  Conference  Cases  374A  and  374B 


Two  plane  strain  flows,  cases  374A  and  374B,  in  1980-81  Stanford 
Conference  were  used  to  further  test  the  zonal  and  K-e  models.  The 
results  of  the  total  and  component  energies  vs.  time  are  shown  in  Figs. 
2.21  and  2.22  for  these  two  cases.  Contrary  to  the  previous  flews,  the 
K-e  model  performs  slightly  better  than  the  zonal  model  for  case  374A 
and  about  equally  well  as  the  zonal  model  for  case  374B.  One  possible 
explanation  is  that  the  initial  conditions  of  experiments  are  aniso¬ 
tropic  which  makes  the  initial  production  rate  not  equal  to  zero  while 
the  zonal  model  starts  with  isotropic  flow  field  and  zero  initial  pro¬ 
duction  rate.  Furthermore,  a  more  careful  study  of  the  experimental 
data  indicates  that  the  Reynolds  stress  structure  in  the  experiments  is 
much  different  from  that  in  Lee  and  Reynolds'  simulation. 

No  model  can  simulate  discrepant  results  for  a  single  physical 
situation.  The  problem  in  this  case  is  that  the  experimental  data  and 
the  full  simulation  are  in  disagreement.  It  is  hard  to  explain  this 
difference.  This  clearly  needs  more  attention.  Until  this  problem  is 
resolved,  these  data  cannot  provide  a  reliable  check  on  model  perfor¬ 
mance.  We  used  the  full  simulations  to  set  the  model  and  are  therefore 
unable  to  predict  the  experimental  data.  The  reverse  could  have  been 
the  case,  i.e.,  if  we  used  the  experimental  data  to  set  the  model,  we 
would  fail  to  predict  the  full  simulation  results. 

2 . 5  Conclusions  of  the  Test 

The  study  of  zonal  models  for  homogeneous  flows,  which  was  meant  to 
test  the  zonal  modeling  concept  and  to  lay  the  ground-work  for  models  of 
more  complex  flows,  has  proven  a  success.  The  results  presented  above 
show  that  zonal  models  perform  far  better  than  the  standard  K-e  model 
in  all  usable  cases.  The  improvement  is  largely  due  to  zonal  models 
being  able  to  avoid  "contamination"  of  the  flow  physics;  of  course,  the 
extra  degrees  of  freedom  which  the  zonal  concept  offers  make  it  much 
easier  to  obtain  accurate  fits  to  the  data.  A  good  model  should  reflect 
the  physics;  the  Reynolds  stress  representation  is  the  tool  for  accomp¬ 
lishing  this.  The  physics  is  not  universal.  This  is  reflected  by  the 
Inability  of  the  Reynolds  stress  function  in  the  standard  K-e  model  to 


fit  the  range  of  flows  considered.  This  is  the  major  drawback  of  the 
standard  K-e  model,  at  least  for  these  flows.  With  the  zonal  approach, 
on  the  other  hand,  we  can  incorporate  more  physics  of  each  zone  into  the 
model  by  constructing  a  Reynolds  stress  function  which  represents  much 
of  the  physics  of  that  zone.  Although  this  requires  more  work  and  a 
careful  study  of  each  flow  zone,  it  can  be  done  once  and  for  all.  The 
gains  obtained  by  whis  work  are  obvious.  As  can  be  seen  in  this  chap¬ 
ter,  the  Reynolds  stress  functions  are  different  for  the  three  homogen¬ 
eous  flows  and  the  ability  to  incorporate  this  is  the  strength  of  the 
zonal  approach. 

Furthermore,  the  zonal  model  has  the  ability  to  adjust  from  one 
zone  to  another  via  a  blending  parameter  governed  by  a  lag  equation. 
The  addition  of  a  lag  equation  provides  one  more  degree  of  freedom. 
Physically,  the  lag  equation  allows  adjustment  of  the  physics  as  the 
flow  changes  its  structure.  This  is  one  of  the  most  attractive  features 
of  zonal  models.  In  this  connection,  it  is  worth  note  that  the  various 
K-e  models,  including  the  standard  model,  as  presented  in  the  1980-81 
Stanford  Conference,  all  use  constants  set  for  equilibrium  situations  in 
homogeneous  flows,  boundary  layers,  and  free  shear  layers.  As  a  result, 
all  these  methods  perform  badly  for  situations  where  flows  are  far  from 
equilibrium  and  in  the  readjusting  regions  that  carry  the  flows  back 
toward  equilibrium.  An  example  is  the  boundary  layer  downstream  of  a 
reattaching  free  shear  layer,  as  occurs  in  the  backward-facing  step. 
Some  additional  degree  of  freedom  In  the  models  is  therefore  apparently 
necessary  for  K-e  type  models  if  they  are  to  provide  good  prediction 
for  nonequilibrium  and  readjusting  flows.  The  zonal  approach  as  used  In 
this  study  supplies  the  necessary  degrees  of  freedom;  thus  the  use  of  a 
"readjustment  model"  is  not  merely  an  artifice  that  allows  fitting,  but 
is  rather  a  necessary  step  in  modeling  the  physics  accurately. 


i 
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In  summary,  the  three  zonal  models  in  this  chapter  were  constructed 
systematically.  They  all  consist  of  a  specific  Reynolds  stress  func¬ 
tion,  a  first-order  lag  equation,  and  two  model  equations  which  have  the 
format  of  the  K-e  model  equations.  The  results  suggest  a  promising 
future  for  the  zonal  modeling  concept. 
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Chapter  III 


A  REVIEW  OF  FREE  SHEAR  FLOW  EXPERIMENTS 


3.1  Preliminary  Remarks 


In  this  chapter,  turbulent  free  shear  flow  measurements  are  re¬ 
viewed.  Numerous  free  shear  flow  experiments  have  been  performed.  How¬ 
ever,  the  agreement  among  experiments  on  the  same  flow  is  rather  poor  in 
some  cases,  particularly  with  respect  to  the  turbulence  data.  Since 
these  data  are  to  be  used  in  constructing  models  and  tuning  model  con¬ 
stants,  it  is  essential  to  justify  the  available  data  and  to  identify 
which  are  reliable. 


There  are  some  reviews  of  f ree-shear-f low  experiments  available  — 
for  example,  Newman  (1966)  and  Harsha  (1971).  However,  they  are  limited 
to  earlier  data  whose  accuracy  is  not  known.  Moreover,  they  deal  mostly 
with  mean-flow  quantities;  very  few  turbulence  quantities  were  examined. 
Rodi  (1972)  made  a  more  thorough  study  of  f ree-shear-f low  experiments 
for  his  modeling  work.  He  reviewed  both  mean  and  turbulence  quantities 
for  a  wide  range  of  flows.  In  a  later  publication,  Rodi  (1975)  provided 
another  detailed  but  similar  survey.  However,  it  has  been  more  than  a 
decade  since  Rodi's  last  review.  Many  more  data  have  recently  become 
available.  Furthermore,  Rodi’s  reviews  contain  no  experiment  that 
covers  both  the  near  and  far  fields  or  a  wide  range  of  the  ratio  of 
velocity  difference  to  free-stream  velocity.  The  details  of  turbulence 
evolution  in  this  whole  range  are  not  clear.  Such  information,  however, 
is  very  important  in  turbulence  modeling;  this  provides  the  motivation 
of  the  present  review. 


Experiments  on  five  basic  free  shear  flows  are  reviewed  in  this 
chapter.  These  flows  are  plane  jets,  plane  wakes  and  mixing  layers, 
axisymmetric  jets,  and  axisymmetric  wakes.  "Pure"  jets  (jets  issuing 
into  stagnant  surroundings)  and  co-flowing  jets  are  considered  sepa¬ 
rately. 


In  addition  to  global  parameters,  such  as  the  spreading  rate  and 
the  centerline  velocity  decay  rate,  we  shall  consider  the  profiles  of 
the  streamwise  mean  velocity,  turbulent  kinetic  energy,  and  turbulent 
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shear  stress  in  this  chapter.  Since  the  turbulent  shear  stress  is  the 
most  important  Reynolds  stress  for  thin  shear  layers,  it  will  receive 
particular  attention.  (Recall  that  it  determines  the  model  constant 
Cp  in  the  K-e  model.) 

In  this  survey,  power  laws  serve  to  correlate  the  data.  The  momen¬ 
tum  integral  equation  provides  a  relationship  between  the  characteristic 
velocity  and  the  width  of  the  layer  and  is  used  as  an  internal  consis¬ 
tency  check.  Moreover,  given  the  streamwise  velocity  profile,  the 
transverse  velocity  and  turbulent  shear  stress  can  be  found  from  the 
continuity  and  momentum  equations.  This  provides  a  check  on  the  shear 
stress  measurements.  There  is  no  easy  way,  however,  to  validate  the 
turbulent  kinetic  energy  data  except  by  comparing  the  profiles  from 
different  experiments.  The  rate  of  dissipation  of  turbulent  kinetic 
energy  is  rarely  measured.  Even  if  measured  and  reported,  the  uncer¬ 
tainties  are  large  and  difficult  to  estimate.  Therefore,  the  dissi¬ 
pation  rate  is  not  examined  here. 


3.2  Plane  Free  Shear  Flows 

3.2.1  Jets 

Many  experiments  have  been  performed  on  plane  jets.  We  shall  con¬ 
sider  some  of  the  well-recognized  ones.  Heskestad  (1965)  and  Gutmark 
and  Wygnanski  (1976)  measured  jets  issuing  into  stagnant  surroundings 
("pure"  jets).  Bradbury  (1965),  in  order  to  reduce  the  measurement 
errors  near  the  edge  of  the  jet,  had  the  jet  exhaust  into  a  slow-moving 
airstream.  Although  this  co-flowing  jet  departs  only  slightly  from  the 
self-preserving  "pure"  jet,  its  spreading  rate  is  considerably  lower 
than  that  of  a  "pure"  jet.  Bradbury  and  Riley  (1967)  carried  out  co¬ 
flowing  jet  experiments  with  various  ratios  of  free  stream  velocity  to 
jet  exhaust  velocity. 

An  empirical  formula  for  the  streamwise  velocity  profile  (Bradbury, 
1965)  fits  the  data  for  all  velocity  ratios  quite  well: 


exp  [~-0.6749n2(l  +  0.0269n4) 


f(n) 


(3.2.1) 


u  -  u. 


(3.2.2) 


Is  the  velocity  excess,  and 


(3.2.3) 


is  the  free  stream  velocity,  Ujm  is  the  maximum  value  of  Uj 
(which  occurs  at  the  centerline),  y  is  the  distance  from  the  center- 
line,  and  $1/2  *8  half-velocity  width  (the  y  value  at  which  Uj  * 
Uim/2).  In  their  "pure"  jet  experiments,  Heskestad  and  Gutmark  and 
Wygnanski  reported  profiles  that  have  a  slightly  less  rapid  approach  to 
the  free  stream  velocity.  However,  the  differences  are  too  small  to  be 
significant.  Furthermore,  Eq.  (3.2.1)  also  represents  the  self¬ 
preserving  profiles  of  axisymmetric  jets  (Bradbury,  1967),  and  wakes 
(Townsend,  1956)  with  only  small  discrepancies.  We  shall  therefore  use 
this  formula  for  all  the  jet  velocity  profiles. 


Co-flowing  jets  have  significantly  different  characteristics  than 
"pure"  jets  although  the  "pure"  jet  is  a  limiting  case  of  the  co-flowing 
jet.  We  shall  examine  co-flowing  jets  separately. 


Pure  Jets 

Jets  issuing  Into  still  surroundings  become  self-similar  after  the 
potential  core  disappears.  The  mean  velocity  reaches  self-similarity 
much  earlier  than  the  turbulence  quantities.  How  fast  this  takes  place 
depends  on  the  nozzle  condition.  In  the  self-similar  region,  the  jet 
width  and  the  centerline  velocity  obey  power  laws.  The  virtual  origin 
differs  from  experiment  to  experiment.  Heskestad  and  Gutmark  and  Wyg¬ 
nanski  suggested  that 


d<$  1/2 

-  0.108  *  0.003 


(3.2.4) 


J  Is  Che  momentum  flux  and  is  Che  invarianc  for  a  pure  jec  by  virtue  of 
Che  momentum  theorem.  Equations  (3.2.4)  and  (3.2.5)  are  also  In  excel¬ 
lent  agreement  with  the  review  of  experiments  by  Rodi  (1972). 

The  turbulent  kinetic  energy  and  shear  stress  profiles  of  various 
experiments  are  compared  in  Figs.  3.1  and  3.2  respectively.  The  agree¬ 
ment  is  within  the  uncertainty  of  the  data.  Among  them,  Bradbury's 
(1965)  data  are  particularly  favored  by  Rodi  (1972).  The  shear  stress 
can  be  obtained  by  solving  continuity  and  momentum  equations  assuming 
the  correctness  of  Eqs.  (3.2.1)  and  (3.2.4);  the  result  is  also  shown  in 
Fig.  3.2.  It  agrees  well  with  Bradbury's  and  Gutmark  and  Wygnanski's 
data.  Heskestad's  data,  on  the  other  hand,  show  a  1 5%  lower  peak  value 
of  the  shear  stress  profile;  these  data  are  not  consistent  with  the 
values  calculated  using  his  own  velocity  profile  (see  Fig.  29  in 
Heskestad,  1965).  He  conclude  that  Heskestad's  data  are  less  reliable 
than  those  of  Bradbury  or  Gutmark  and  Uygnanski. 

Co-Flowing  Jets 

A  jet  issuing  into  a  uniformly  moving  stream  cannot  be  self-similar 
because  U^/Uj  varies  with  x  (Townsend,  1956);  therefore  we  are 
interested  in  the  streamwise  flow  development.  To  our  knowledge,  only 
Bradbury  and  Riley  (1967)  reported  co-flowing  jet  measurements  cover¬ 
ing  a  wide  range  of  velocity  ratios.  They  showed,  using  dimensional 
arguments,  that  the  co-flowing  jet  approaches  the  behavior  of  a  self- 
preserving  "pure”  jet  as  (x-xQ)/8  -►  0  and  a  self-preserving  wake  as 
(x-xQ)/0  ♦  ».  They  verified  these  results  experimentally.  Here,  xQ 
is  the  apparent  origin  of  the  flow  and  is  the  only  influence  the  nozzle 
conditions  have  on  the  flow  and  9  is  the  momentum  thickness  of  the  jet 
defined  as 


The  limiting  behaviors  were  also  suggested  by  Townsend  (1956). 


Bradbury  and  Riley  reported  that  the  data  can  be  correlated  in  two 
limiting  cases  as  follows.  In  the  limit  (x-xq)/0  ♦  0, 
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(3.2.8a) 


and 


-  0. 16 


In  the  limit  (x-xq)/9  +  », 


and 


-  0.41 
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(3.2.8b) 


(3.2.9a) 


(3.2.9b) 


Note  that  Eqs.  (3.2.8a)  and  (3.2.8b)  reduce  to  Eqs.  (3.2.4)  and  (3.2.5) 
as  Uw  +  0;  this  shows  that,  in  the  early  stage,  a  co-flowing  jet 
behaves  similarly  to  a  "pure"  jet.  In  addition,  Bradbury  and  Riley 
provided  the  width  and  centerline  velocity  excess  versus  downstream  dis¬ 
tance  in  a  dimensionless  form  that  collapses  all  the  co-flowing  jet  data 
except  for  the  initial  developing  stage,  cf..  Figs.  3.3  and  3.4. 


With  the  mean  velocity  profile  represented  by  Eq.  (3.2.1),  the  mo¬ 
mentum  integral  equation,  Eq.  (3.2.7),  becomes 
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Equations  (3.2.10)  and  (3,2.11)  provide  a  relationship  between  the  jet 
width  and  the  centerline  velocity  excess.  In  the  limit  (x-xQ)/0  *  0 
(self-preserving  "pure"  jet),  the  first  terra  on  the  left-hand-side  of 
Eq.  (3.2.10)  dominates.  On  the  other  hand,  in  the  limit  (x-xQ)/0  ♦  « 
(self-preserving  wake),  the  second  term  dominates.  Therefore,  Eq. 
(3,2. 10)  becomes 
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respectively,  in  the  two  limits, 


Equations  (3.2.8)  and  (3.2.9)  satisfy  Eqs.  (3.2.12)  and  (3.2.13), 
respectively,  within  2X,  providing  a  consistency  check.  Furthermore, 
co-flowing  jet  data  given  by  Bradbury  and  Riley  also  agree  with 
Eq.  (3.2.10)  within  2X. 


Neither  turbulent  kinetic  energy  nor  shear  9tress  profiles  were 


reported  by  Bradbury  and  Riley. 


However,  the  eddy  Reynolds  number 


defined  by 
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(3.2. 14) 


was  tabulated.  The  only  conclusion  that  can  be  drawn  from  the  eddy 


Reynolds  number  data  is  that,  when  going  downstream  in  a  co-flowing  jet, 

2 

the  shear  stress  level  (normalized  by  U  )  increases  to  an  asymptotic 

itn 


upper  limit  equivalent  to  far-field  wake  flows.  This  Indicates  that  a 
far-field  co-flowing  jet  behaves  similarly  to  a  far-field  wake.  These 
results  will  help  us  clarify  the  controlling  physical  parameters  in  free 


shear  flows. 


3.2.2  Wakes 


Plane  wake  experiments  include  those  behind  flat  plates,  those 
behind  circular  cylinders,  and  some  intermediate  cases.  Although  the 
near  fields  of  wake  flows  depend  heavily  on  the  initial  conditions, 
i.e.,  the  shape  of  the  body  that  generates  the  wake  and  whether  or  not 
there  is  vortex  shedding,  the  far  fields  display  a  universal  asymptotic 
self-similar  state.  Townsend  (1949)  measured  a  wake  behind  a  circular 
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cylinder  and  claimed  that  self-similarity  is  achieved  500  diameters 
downstream  of  the  cylinder.  The  centerline  velocity  deficit  at  the  last 
station  was  about  32  of  the  free  stream  velocity.  A  complete  set  of 
turbulence  measurements  was  reported.  However,  Townsend's  data  were 
taken  only  in  the  far-field  asymptotic  self-similar  region.  They  do 
not,  therefore,  provide  any  information  on  the  evolution  of  the  flow 
into  the  asymptotic  state.  Chevray  and  Kovasznay  (1969),  Andreopoulos 
and  Bradshaw  (1980),  and  Ramaprian,  Patel  and  Sastry  (1982)  investigated 
the  wake  of  a  flat  plate.  Their  measurements  were  made  in  the  near 
field;  the  centerline  velocity  deficit  at  the  last  measurement  station 
is  at  least  20%  of  the  free  stream  velocity.  An  analysis  of  the  momen¬ 
tum  integral  equation  similar  to  Eq.  (3.2.10)  and  examining  the  size  of 
the  terras  show  that  these  flows  are  not  near  the  far  wake  asymptote. 

More  recently,  Pot  (1979)  studied  the  flat  plate  wake  flow  from  the 
trailing  edge  of  the  plate  to  a  location  sufficiently  far  downstream 
that  the  far  wake  asymptote  is  clearly  reached.  He  reported  all  the 
turbulence  data  necessary  for  the  present  study.  Pot's  data  set  and 
that  of  Townsend  will  be  examined  here. 

The  far-field  velocity  profiles  of  Townsend  and  Pot  are  well  ap¬ 
proximated  by  Eq.  (3.2.1)  when  Uj  is  taken  to  be  the  velocity 
deficit  Ua  -  U.  Velocity  profiles  at  the  last  few  stations  of  Chevray 
and  Kovasznay,  Andreopoulos  and  Bradshaw,  and  Ramaprian,  Patel  and 
Sastry  all  agree  with  this  empirical  formula  although  their  experiments 
were  not  carried  far  enough  downstream  to  fully  reach  the  asymptotic 
state.  This  indicates  that  the  mean  flow  quantities  settle  down  to  a 
self-similar  stage  much  earlier  than  the  turbulence  quantities,  as  has 
been  found  in  a  wide  variety  of  other  flows. 


Pot's  far-field  data  suggest 
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These  equations  agree  with  (3.2.9a)  and  (3.2.9b)  which  describe  the 
behavior  of  the  far-field  co-flowing  jet.  This  again  confirms  that  far- 
field  wakes  and  far-field  co-flowing  jets  have  common  characteristics. 
A  spreading  parameter  defined  by 


U  dfi, 


(3.2.17) 


is  0.103  according  to  Eqs.  (3.2.15)  and  (3.2.16),  compared  to  a  value  of 
0.098  reported  by  Rodi  (1972)  from  Townsend's  data.  The  5%  difference 
is  within  the  uncertainty  of  the  data. 


The  counterpart  of  Eq.  (3.2.10)  for  wake  flow  is 
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Pot's  wake  width  and  centerline  velocity  deficit  data  were  checked 
against  Eq.  (3.2.18).  The  agreement  is  excellent  (within  3%)  except  at 
the  first  two  stations  where  the  velocity  profiles  are  not  close  to  the 
self-similar  form. 

It  is  clear  from  Pot's  data  that  both  the  turbulent  kinetic  energy 
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and  shear  stress  (both  normalized  by  )  increase  with  downstream 
distance  to  the  asymptotic  far  wake  values.  Figure  3.5  shows  the  maxi¬ 
mum  shear  stress  and  maximum  kinetic  energy  in  the  layer  as  functions  of 
downstream  distance.  The  far-field  shear  stress  level  is  between  0.05 
and  0.053,  which  agrees  well  with  Townsend's  value  (0.051).  Moreover, 
using  0.103  for  the  spreading  parameter  defined  by  Eq.  (3.2.17)  and  the 
velocity  profile,  Eq.  (3.2.1),  the  shear  stress  can  be  derived  from  the 
momentum  equation.  The  resulting  maximum  shear  stress  is  about  0.053 
which  agrees  with  the  data  from  both  sources.  The  far-field  kinetic 
energy  level  is  about  0.18  in  Pot's  data  and  less  than  0.14  in  Town¬ 
send's  data.  It  is  hard  to  justify  the  discrepancy  between  these  two 
data  sets.  However,  we  tend  to  believe  Pot's  data,  for  the  following 
reasons.  Firstly,  Pot's  experiment  is  much  more  recent,  and  his  equip¬ 
ment  is  more  modern.  Since  much  has  been  learned  about  hot-wire  error 
control  during  the  intervening  period,  he  should  get  more  reliable 
results.  Secondly,  all  known  plane  shear  flows  have  |u'v'|/k  =  0.3  in 


the  equilibrium  region.  In  Townsend's  data  this  ratio  become  0.5,  while 

for  Pot's  data  it  remains  at  about  0.3.  Thirdly,  similar  experiments 

performed  by  Uberoi  and  Freymuth  (1969)  and  Thomas  (1973)  suggested  that 

2 

Townsend's  values  of  u'  are  too  snail. 

3.2.3  Mixing  Layers 

The  turbulent  plane  nixing  layer  is  one  of  the  simplest  shear 
flows.  However,  complete  understanding  of  the  structure  of  this  flow  Is 
still  far  from  complete  although  the  general  features  of  the  flow  were 
established  by  Liepmann  and  Laufer  (1947).  They  made  hot  wire  mea¬ 
surements  in  a  zero-velocity-ratio  mixing  layer.  Their  mean  velocity 
and  turbulence  measurements  indicated  a  very  rapid  approach  to  self¬ 
similarity.  This  was  regarded  as  the  most  reliable  and  extensive  source 
of  mixing  layer  data  for  more  than  two  decades.  Wygnanski  and  Fiedler 
(1970),  Patel  (1973),  and  Champagne,  Pao  and  Wygnanski  (1976)  studied 
the  same  flow  to  obtain  higher-order  statistical  characteristics  of  the 
turbulence.  The  agreement  among  these  experiments  is  rather  poor  as 
will  be  noted  later. 

Two-stream  mixing  layers  with  finite  velocity  ratios  have  not  re¬ 
ceived  as  much  attention.  Sabin  (1965)  measured  mean  velocity  profiles 
and  correlated  the  spreading  parameter  for  various  velocity  ratios  with 
a  simple  function.  Miles  and  Shih  (1968)  investigated  the  flow  with  a 
wide  range  of  velocity  ratios.  However,  only  the  spreading  parameter  as 
a  function  of  velocity  ratio  was  presented.  Spencer  and  Jones  (1971) 
and  Yule  (1971)  each  measured  the  layers  at  two  different  velocity 
ratios.  Both  reported  detailed  velocity  profiles  and  turbulence  quan¬ 
tities.  The  agreement  is  again  poor. 

Mixing  layers  are  known  to  be  very  sensitive  to  the  initial  condi¬ 
tions.  For  example,  placing  a  trip  wire  upstream  on  the  splitter  plate 
can  alter  the  spreading  rate  by  more  than  20%  (Liepmann  and  Laufer,  1947 
and  Wygnanski  and  Fiedler,  1970).  They  are  also  sensitive  to  the  bound¬ 
ary  conditions  imposed.  In  most  cases,  mixing  layer  experiments  are 
conducted  in  wind  tunnels  with  walls  both  above  and  below  the  flow. 
This  affects  the  secondary  flow  speed,  the  entrainment  into  the  mixing 
layer  and,  consequently,  the  turbulence  structure  and  the  spreading 


35 


% 


:«sa 


k 


jSj 

v? 

$ 


9 


i 


rate.  Variations  In  Initial  and  boundary  conditions  are  probably 
responsible  for  the  disagreement  among  experimental  data. 


Figures  3.6  compares  velocity  profiles  for  mixing  layers  with  zero 
and  nonzero  velocity  ratios.  In  this  figure,  y  is  normalized  by  a 


characteristic  flow  width  6  -  y»  Q  -  y„  . .  Here,  y  is  defined  as 

U#  7  U •  1  X 


the  location  where  U  -  U_  ♦  x(U_  -  U_  ).  There  Is  considerable  scat- 

2  1  2 

ter  of  the  data  near  the  edges,  especially  at  the  zero  velocity  edge  in 


Fig.  3.6.  According  to  Rod!  (  1975),  however,  this  is  the  best  way  to 
correlate  the  data.  If  the  data  are  plotted  versus  y/x,  as  is  often 
done  in  the  literature,  the  agreement  is  worse  because  the  spreading 
rate  dd/dx  varies  considerably.  The  scatter  Indicates  the  diflculties 
in  measuring  accurately  in  these  regions  and  perhaps  reflects  the  sensi¬ 
tivity  of  the  layer  to  the  nearness  of  the  boundary  walls. 


Figure  3.7  shows  spreading  rates,  dL/dx,  as  a  function  of  the 


velocity  ratio  (R  -  U  /U  ).  Note  that  the  layer  width  L  in  this 

“2  “1 

figure  is  defined  as 


y  _ -  y _ 

✓0.9  /0.1 


Some  of  the  data  were  taken  from  Rodi  (  1975).  The  line  in  the  figure 
represents  the  spreading  rate  given  by  the  Evaluation  Committee  of  1980- 
81  Stanford  Conference.  When  the  extreme  values  are  discarded,  the  line 
seems  a  reasonable  average  representation  for  the  spreading  rate, 
although  it  is  slightly  lower  than  the  average  at  high  velocity  ratios. 


The  profiles  of  turbulent  kinetic  energy  measured  in  various  exper¬ 
iments  are  shown  in  Fig.  3.8(a)  and  (b)  for  zero  and  two  finite  velocity 
ratios  (R  “  0.3  and  0.61).  The  agreement  is  poor  except  for  the 
shape  of  the  profiles. 


The  turbulent  shear  stress  profiles  are  shown  in  Fig.  3.9(a)  and 
(b)  for  zero  and  two  finite  velocity  ratios  (R  -  0.3  and  0.61), 
respectively.  Although  the  agreement  is  not  good,  it  appears  that  the 
nondimensional  shear  stress  increases  as  the  velocity  ratio  is  in¬ 


creased. 


There  is  no  way  to  justify  the  disagreement  because  of  the  dif¬ 
ferent  experimental  setups  in  different  experiments  and  the  nature  of 


%-v 


WWW* 


V  aVj  4  ' 


sensitivity  of  the  flow  as  explained  above.  Until  this  issue  is  com¬ 
pletely  resolved,  the  only  quantity  that  can  be  used  as  a  guide  for 
modeling  seems  to  be  the  average  spreading  rate  given  by  1980-81 
Stanford  Conference. 

3.2.4  Shear  Stress-Kinetic  Energy  Relationship 

The  value  of  |u'v'|/K  was  examined  for  all  plane  shear  flows. 
For  jets  and  wakes,  this  quantity  increases  linearly  from  zero  at  the 
centerline  to  a  maximum  of  about  0.3  and  then  drops  back  to  zero. 
Townsend's  (1949)  wake  flow  is  an  exception;  it  shows  a  maximum  value 
much  higher  than  0.3.  However,  it  is  believed  that  his  measurements  of 
kinetic  energy  are  too  low  as  explained  above.  For  mixing  layers,  all 
the  available  data  show  that  |u'v'|/K  has  a  value  of  about  0.3  in  the 
center  and  decreases  to  zero  at  both  edges. 

3. 3  Axlsymmetrlc  Free  Shear  Flows 

3.3.1  Jets 

The  mean  flow  field  of  axisymmetric  jets  have  been  studied  by  nu¬ 
merous  researchers  (for  example,  Hinze,  1959),  but  the  turbulence  field 
by  only  a  few.  Most  of  the  literature  on  this  flow  appeared  more  than 
three  decades  ago  when  measuring  techniques  were  not  very  advanced. 
More  recent  measurements  were  made  by  Wygnanski  and  Fiedler  (1969)  and 
Rodi  (  1972)  for  self-preserving  pure  jets  and  by  Maczynski  (1962), 
Reichardt  (1965),  and  Antonia  and  Bilger  (1973)  for  co-flowing  jets.  It 
is  again  easier  to  discuss  them  separately. 

Pure  Jets 

It  was  suggested  by  Bradbury  (  1967)  that  Eq .  (3.2.1)  can  also  ap¬ 
proximate  the  velocity  profiles  of  axisymmetric  jets.  However,  the  pure 
jet  velocity  profiles  of  both  Wygnanski  and  Fiedler  (1969)  and  Rodi 
(1972)  approach  zero  more  slowly  than  the  formula  suggests.  A  similar 
trend  was  found  In  plane  jet  cases.  However,  the  uncertainty  of  the 
measurements  is  large  near  the  edge  (especially  when  the  edge  has  zero 
velocity).  Consequently,  this  empirical  formula  suffices  as  a  first 
order  approximation. 
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The  spreading  rates  of  the  pure  jets  from  both  Wygnanski  and 
Fiedler  (1969)  and  Rodi  (1972)  are 


0.087  ±  0.002 


(3.3.1) 


This  value  agrees  with  the  spreading  rate  reported  by  Newman  (1967), 
d6i/2/dx  ■  0.086,  an  average  over  earlier  experiments. 


The  centerline  velocity  decay  rate,  according  to  Wygnanski  and 
Fiedler,  is 


d(VUlm) 


d(x/D) 


-  0.2 


in  self-preserving  region,  where  Uj  is  the  jet  nozzle  exhaust  veloc¬ 
ity.  This  converts  into 


-  0.131 


(3.3.2) 


where 
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is  the  momentum  flux;  it  is  the  invariant  for  the  jet. 


The  turbulent  kinetic  energy  and  shear  stress  profiles  of  Wygnanski 
and  Fiedler  and  Rodi  are  presented  in  Figs.  3.10  and  3.11,  respectively. 
The  agreement  is  reasonable.  The  maximum  non-dimensional  shear  stress 
of  Wygnanski  and  Fiedler  is  about  10%  less  than  Rodi's  value.  The  lat¬ 
ter  is  closer  to  the  value  (0.0184)  calculated  from  the  continuity  and 


momentum  equations. 


Figure  3.12  shows  the  ratio  Ju'v'l/K  versus  the  transverse 
distance  for  axisymmetric  and  plane  pure  jets.  The  maximum  value  of 
|u'v'|/K  for  the  axisymmetric  cases  is  about  0.23  vs.  0.3  for  the 
plane  cases.  We  speculate  that  this  is  due  to  lateral  vortex  stretching 
which  appears  only  in  axisymmetric  flows. 
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Co-Flowing  Jets 

The  velocity  profiles  measured  by  Antonia  and  Bilger  (1973)  agree 
with  Eq.  (3.2.1)  satisfactorily,  as  reported  by  Bradbury  (1967).  Rodi 
(1975)  also  claims  that  the  shape  of  velocity  profile  is  approximately 
the  same  for  both  plane  and  axisymraetric  jets,  and  that  it  does  not  vary 
significantly  with  either  velocity  ratio  or  downstream  distance. 


Before  discussing  the  spreading  and  centerline  velocity,  it  is  wise 
to  examine  the  constraint  connecting  them.  With  the  momentum  thickness 
for  axisymmetric  jets  defined  by 


and  the  help  of  Eq.  (3.2.1),  we  can  write 


(3.3.4) 
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In  the  case  of  a  pure  jet,  Eq.  (3.3.5)  reduces  to 
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The  pure  jet  described  by  Eqs.  (3.3.1)  and  (3.3.2)  satisfies  Eq.  (3.3.7) 
within  4%. 

The  spreading  and  velocity  decay  rates  in  Antonia  and  Bilger's 
experiment  show  large  uncertainty  and  scatter.  Furthermore,  the  data  do 
not  display  internal  consistency.  Rodi  (1975)  managed  to  collapse  the 
jet  width,  Fig.  3.13,  and  centerline  velocity  data.  Fig.  3.14,  of  three 
experiments  for  a  wide  range  of  velocity  ratios.  At  larger  values  of 
x/0,  however,  the  measured  jet  widths  and  centerline  velocities  do  not 
satisfy  Eq.  (3.3.5).  Since  the  measurements  of  centerline  velocity  are 
probably  more  accurate  than  those  of  the  jet  width,  a  jet  width  computed 
from  Eq .  (3.3.5)  with  the  measured  centerline  velocity  given  is  also 
shown  in  Fig.  3.13. 
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No  measurements  of  turbulent  kinetic  energy  or  shear  stress  were 

taken  by  Maczynski  (1962)  or  Reichardt  (1965).  Antonia  and  Bilger's 

(1973)  data  show  that  both  the  scaled  kinetic  energy  and  shear  stress 

Increase  with  downstream  distance  In  their  measuring  range.  The  largest 
-  2 

value  of  |u'v'|/U  they  recorded  is  about  three  times  that  found  in 

lm 

pure  jets.  Although  their  data  must  be  considered  unreliable  because  of 
the  scatter,  the  increasing  trend  of  both  the  non-dimensional  kinetic 
energy  and  shear  stress  in  co-flowing  jets  seems  to  be  established. 
This  trend  was  also  found  in  the  plane  cases.  Indeed,  the  shear  stress 
obtained  from  the  continuity  and  momentum  equations  and  the  assumed 
velocity  profile  confirms  this. 

With  the  assumption  that  the  velocity  fluctuations  in  the  radial 
direction  are  about  the  same  in  magnitude  as  those  in  the  streamwise 
direction,  Antonia  and  Bilger's  turbulence  data  show  a  maximum  value  of 
|u'v'|/K  about  0.3.  Since  their  data  were  taken  some  distance  from 
the  nozzle,  where  the  spreading  rate  is  considerably  lower  than  the  pure 
jet  spreading  rate,  this  leads  us  to  believe  that  the  effect  of  vortex 
stretching  is  significant  only  when  the  jet  spreads  rapidly.  We  shall 
use  this  fact  in  generalizing  our  parameterization  of  free  shear  flows. 

3.3.2  Wakes 

Axisymmetric  wakes  are  known  to  be  very  sensitive  to  the  shape  of 
the  wake-generating  body  (Rodi,  1975).  Unless  the  body  is  very  slender 
it  usually  creates  vortices  in  Its  wake.  Different  body  shapes  generate 
vortices  of  different  strengths.  This  in  turn  provides  different  ini¬ 
tial  conditions  for  the  various  wakes  and  it  takes  a  very  long  distance 
for  the  flow  to  "forget"  the  effects  of  the  initial  conditions.  Al¬ 
though  asymptotic  far-field  behavior  exists,  it  usually  occurs  very  far 
downstream  from  the  body. 

Many  axisymmetric  wakes  have  been  studied  (e.g.,  Carraody,  1964; 
Chevray,  1968;  and  Uberoi  and  Freymuth,  1970).  Most  of  them  are  gen¬ 
erated  by  blunt  bodies,  e.g.,  cones,  spheres  or  disks.  Unfortunately, 
they  were  not  studied  for  a  long  downstream  distance  because  the  re¬ 
searchers  were  Interested  mainly  in  the  vortex  shedding  phenomenon. 
Only  Chevray  (1968)  measured  a  wake  flow  that  was  extended  close  to  the 
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self-preserving  asymptote.  His  wake-generating  body  was  a  relatively 
slender  6:1  spheroid  and  the  effect  of  vortex  shedding  was  reduced.  A 
complete  set  of  mean  flow  and  turbulence  quantities  was  reported. 


Chavray's  self-similar  velocity  profiles  are  well  approximated  by 
Eq.  (3.2.1).  The  centerline  velocity  deficit  and  the  wake  width  both 
approach  asymptotic  power  laws.  Moreover,  the  mean  flow  quantities 
satisfy  (to  within  3%)  the  constraint 


(3.3.8) 


which  can  be  derived  in  manner  of  Eq .  (3.3.5).  According  to  Rodi 
(1975),  the  far-wake  spreading  parameter  defined  by  Eq.  (3.2.17)  for 
this  flow  is  about  0.105,  the  same  as  that  of  Reichardt's  (1965)  asymp¬ 
totic  co-flowing  jet.  This  confirms  that,  as  in  plane  cases,  an  axi- 
symmetric  co-flowing  jet  behaves  similarly  to  an  axisymmetric  wake  in 
the  far  field. 

The  turbulence  field  of  Chevray's  wake  exhibits  a  trend  similar  to 
that  seen  in  co-flowing  jets.  The  nondimensional  kinetic  energy  and 
shear  stress  increase  to  asymptotic  values  as  the  flow  develops.  Fig. 
3.15  shows  the  maximum  values  of  kinetic  energy  and  shear  stress  in  the 
layer  as  functions  of  downstream  distance. 


3 . A  Conclusions 

Experiments  on  five  basic  free  shear  flows  were  reviewed  in  this 
chapter,  including  both  near  and  far  fields,  the  approach  to  asymptotic 
behavior,  and  co-flowing  cases.  Whenever  possible,  the  data  were 
checked  by  analytical  relationships  for  consistency.  Good  data  were 
identified  for  each  flow.  For  those  cases  where  scatter  and  uncertainty 
were  large,  trends  were  established;  this  will  help  us,  at  least  quali¬ 
tatively,  to  understand  the  physics.  It  appears  that  the  measurements 
become  difficult  when  the  velocity  or  velocity  difference  is  small,  for 
example,  near  the  edges  of  shear  layers.  This  is  probably  due  to  the 
fact  that  the  uncertainty  of  the  data  is  larger  than  the  measured  val¬ 
ues  . 
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In  later  chapters,  we  shall  use  the  data  reviewed  in  this  chpater 
to  construct  models  and  tune  the  model  constants.  The  conclusions  about 
limiting  forms  and  their  relations  to  each  other  through  nondimensional 
correlations  established  in  this  chapter  will  allow  us  to  create  a 
single  model  for  free  shear  flows  containing  two  nondimensional  param¬ 
eters  and  a  simple,  uniform  method  for  treating  readjustment  regions. 


8ft 
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Chapter  IV 
NUMERICAL  METHOD 


4. 1  Introduction 

As  a  result  of  the  thin  shear-layer  approximation,  the  governing 
equations  for  both  the  mean  and  turbulence  fields  for  free  shear  flows 
are  parabolic.  They  are: 
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where  vt  =  C^K  /e  is  the  eddy  viscosity.  In  these  equations,  1=0 
Is  for  plane  flows  and  i  =  1  for  axisymmetric  flows.  Note  that  the 
K-e  model  equations  are  used  in  Eqs.  (4.1.3)  and  (4.1.4)  for  the  turbu¬ 
lence  field.  The  standard  method  of  solving  parabolic  equations  is  to 
march  the  solution  downstream  in  the  predominant  flow  direction.  This 
Is  because  Information  propagates  in  the  flow  direction  and  what  happens 
at  any  location  has  no  effect  upstream.  This  behavior  is  a  consequence 
of  the  characteristics  of  these  equations.  In  the  flow  direction,  the 
domain  is  open  and  only  an  initial  upstream  condition  is  needed.  Two 
boundary  conditions  are  required  in  the  transverse  direction. 


The  Keller  Box  method  (see,  for  example,  Cebeci  and  Bradshaw,  1977) 
is  a  popular  method  for  solving  parabolic  equations.  This  method 
approximates  all  variables  by  their  values  at  points  on  a  rectangular 
grid  and  replaces  the  differential  equations  by  finite  difference 
equations  that  are  averaged  over  each  cell  (or  box).  Figure  4.1(a) 
shows  a  typical  grid  box  used  in  this  method.  However,  experience  shows 
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that  this  method  tends  to  give  oscillatory  solutions  in  the  marching 
direction  if  the  initial  condition  is  not  perfect.  If  the  initial  con¬ 
dition  is  too  high,  the  solution  at  the  first  step  tends  to  be  low  so 
that  the  average  at  the  center  is  correct;  similarly,  the  solution  at 
the  second  step  is  high.  The  oscillation  continues  downstream  although 
its  magnitude  reduces.  Another  drawback,  of  this  method  is  its  require¬ 
ment  of  a  V  velocity  component  initial  condition.  This  component  is 
rarely  measured  in  experiments  and  has  to  be  estimated  from  the  continu¬ 
ity  equation. 

This  provides  the  motivation  for  the  development  of  a  better  numer¬ 
ical  scheme  for  solving  the  governing  equations  for  free  shear  flows. 
In  this  chapter,  a  modified  Keller  Box  method  that  uses  a  staggered  grid 
system  is  developed.  A  transformation  of  the  governing  equations  equiv¬ 
alent  to  using  a  grid  expanding  with  the  growth  of  the  layer  is  first 
performed.  An  averaging  technique  at  the  first  step  of  the  solution 
procedure  eliminates  oscillations.  Variable  step  sizes  determined  by 
the  growth  of  the  layer  are  used.  Finally,  the  method  is  tested  for 
some  standard  cases.  The  method  of  generation  of  initial  conditions  is 
described. 

4 . 2  Description  of  the  Numerical  Scheme 

4.2.1  Transformation  of  the  Governing  Equations 

To  account  for  the  growth  of  the  layer,  the  governing  equations  are 
transformed  from  the  (x,y)  physical  coordinates  to  (C,n)  variables 
where 


£  -  x 


(4.2.  1  ) 


and 


n 


(4.2.2) 


Here,  E(x)  is  a  constant  proportional  to  the  thickness  of  the  shear 
layer.  Thus,  In  the  transformed  coordinates,  the  layer  thickness  is 
constant  and  the  same  number  of  grid  points  can  be  used  at  each  down¬ 
stream  location. 
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After  transformation,  the  governing  equations,  Eqs.  (4.1.1)  to 
(4. 1.4)  become 
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Before  solving  Eqs.  (4.2.3) 
fied.  They  are  determined, 
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to  (4.2.6),  E  and  dE/dx  must  be  speci- 
in  a  discrete  manner,  as  follows: 
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where  the  subscript  represents  the  step  number  in  the  streamwlse  direc¬ 
tion.  dE  dx  can  be  calculated  once  C  Is  determined.  The  half-velocity 
width  (5  ,  /,)n  Is  easily  determined  from  the  converged  solution  at 


4  5 


■»  .% 


.-  /7.V.V.V.  a 


r ,  -r 


4.2.2  Grid  System 


.y 

v 

v 


J 


To  derive  the  finite  difference  equations  for  a  box,  it  is  first 
necessary  to  define  the  locations  at  which  the  variables  are  defined. 
As  shown  in  Fig.  4.1(b),  U,  K,  and  e  are  given  at  the  centers  of 
left  and  right  boundaries  and  V  at  the  centers  of  top  and  bottom  bound¬ 
aries  of  a  box.  This  arrangement  makes  the  streamwise  derivatives  of 
U,  K  and  e,  and  normal  derivative  of  V  much  easier  to  deal  with. 
Furthermore,  no  V  velocity  component  need  be  provided  at  the  initial 
station. 


The  outermost  location  of  U,  K  and  e  is  outside  the  last  box. 
Therefore  the  zero  boundary  conditions  of  both  K  and  e  on  the  edge 
do  not  cause  any  difficulty  in  the  calculation  of  quantities  such  as 
K^/e  or  e^/K  in  the  last  box. 


The  spacing  between  grid  points  in  the  normal  direction  can  be 
either  constant  or  variable.  The  choice  is  made  by  the  selection  of  the 
initial  grid.  Details  will  be  given  in  next  section.  The  spacing  in 
the  streamwise  direction  is  variable,  and  the  increase  in  step  size  is 
determined  by  the  rate  at  which  the  layer  grows: 
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where  Ax  =  x  .  ,  -  x  . 

n  n  +  1  n 


4.2.3  Formulation  of  the  Finite  Difference  Equations 

Central  differencing  is  used  everywhere  except  where  the  streamwise 
velocity  is  small  relative  to  the  normal  velocity.  This  occurs  near  the 
edge  of  pure  |ets.  For  those  cases,  the  flow  is  principally  in  the 
normal  direction  near  the  edge  and  central  differencing  causes  oscilla¬ 
tion  when  the  cell  Reynolds  number  (based  on  V  and  Ay)  is  larger  than 
two  (see  Patankar,  1980).  To  remedy  this,  upwind  differences  must  be 
used  for  the  convection  terms  in  the  normal  direction  in  this  region 
(Patankar,  19801.  A  second-order  upwind  scheme  is  used  near  the  edge. 

The  fomulation  of  the  flnit<*  difference  equations  is  similar  to 
that  of  the  Keller  Box  method  and  can  he  found  in  Cebec  i  and  Bradshaw 
! I  9 7 7  ) .  "he  details  are  presented  in  Appendix  A. 
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4.2.4  Averaging  at  the  First  Ste 


To  avoid  oscillation  of  the  solution  in  the  streamwise  direction, 
an  averaging  technique  is  used  at  the  first  step.  The  initial  condi¬ 
tions  are  generally  not  correct.  The  finite  difference  equations 
represent  an  average  of  the  governing  equations  and  tend  to  produce 
oscillations  with  an  amplitude  equal  to  the  initial  error.  Averaging  of 
the  solutions  at  the  initial  and  first  stations  eliminates  most  of  the 
error  and  provides  an  Improved  initial  condition.  The  solution  proce¬ 
dure  is  then  restarted  from  the  center  location  and  a  half  step  is  taken 
so  that  an  improved  solution  is  obtained  at  the  first  station.  This 
averaging  process  is  repeated  three  times.  Each  time  the  location  of 
averaged  solution  is  moved  forward  by  half  the  reduced  step  size. 

After  this  averaging,  the  effect  of  "bad"  initial  conditions  on  the 
solution  is  greatly  reduced.  From  the  second  step  on,  the  normal  solu¬ 
tion  procedure  is  used. 


4.2.5  Calculation  of  Eddy  Viscosit) 


The  eddy  viscosity  is  required  before  the  system  of  equations  can 
be  solved.  In  order  to  calculate  the  eddy  viscosity,  K  and  e  are 
needed.  Since  the  solution  procedure  is  iterative,  K.  and  e  are 
lagged  by  one  iteration  in  the  calculation  of  eddy  viscosity.  No  other 
variables  are  lagged  in  the  computation;  mean  and  turbulence  quantities 
are  computed  simultaneously. 

Near  the  edge  of  the  layer,  the  profile  of  eddy  viscosity  is 
smoothed  when  necessary.  Bad  profiles  are  sometimes  caused  by  small 
negative  values  of  K  or  e.  Because  these  quantities  are  very  small 
in  this  region,  they  may  produce  large  errors  in  the  eddy  viscosity. 
However,  eddy  viscosity  snould  approach  zero  near  the  edge  of  the 
laver.  Therefore,  whenever  a  bad  value  of  the  eddy  viscosity  occurs,  a 
/alue  linearly  interpolated  between  the  previous  good  point  and  zero  at 
the  edge  replaces  the  bad  value.  This  smoothing  process  prevents  the 
had  solution  from  penerratlng  into  the  layer.  It  is  never  necessary  to 
use  this  procedure  for  more  than  three  or  four  points  near  the  edge 
except  ‘‘or  far  wake  flows. 


The  solution  method  Is  Iterative.  The  convergence  criterion 
demands  that  the  increments  of  U,  K,  and  e  at  the  point  nearest  the 
centerline  between  iterations  be  less  than  10-^  times  the  values  of  the 
quantities  at  the  previous  iteration.  Typically,  the  solution  converges 
within  five  iterations. 

4.2.7  Solution  Procedures 

Provided  the  solutions  at  step  n  have  been  obtained,  the  solution 
procedure  can  be  outlined  as  follows: 

1.  In  physical  coordinates,  calculate  &\/2  at  steP  n.  Calculate 
En+j  by  Eq.  (4.2.8).  Calculate  dE/dx. 

2.  Use  the  solutions  at  step  n  as  initial  guesses  for  all  variables  at 
step  n+1. 

3.  Solve  a  block  tridiagonal  system  of  linear  equations  that  results 
from  the  Newton-Raphson  linearization  of  the  finite  difference 
equations . 

4.  Update  all  variables  by  adding  the  increments  obtained  from  step  3 
to  the  old  values. 

5.  Update  the  eddy  viscosity. 

6.  If  the  Increments  satisfy  the  convergence  criteria,  stop  the  Itera¬ 
tion  and  go  to  the  next  station.  If  not,  repeat  steps  3  to  5. 


4.3  Input  Data 


4.3.1  Initial  Conditions 

Initial  profiles  of  U,  K  and  e  are  required  to  start  the  com¬ 
putation.  They  have  to  be  generated  from  experimental  data.  The 

streamwise  velocity  is  measured  in  all  experiments.  For  turbulent 

2  2 

kinetic  energy  components,  u'  and  v'  are  measured  in  most  exper- 
2 

lments.  w'  ,  on  the  other  hand,  is  not  often  measured  and  is  usually 

F  — 2 

approximated  by  the  average  of  u’  and  v*  The  dissipation  rate  of 

■ 

turbulent  kinetic  energy  is  rarely  measured;  if  it  is,  the  uncertainty 
is  very  high.  The  turbulent  shear  stress  is  always  recorded.  One  way 
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to  calculate  the  dissipation  rate  is  via  the  relationship  for  the  shear 


stress  used  in  the  K-e  model: 


U  (-u’v') 


(4.3.1) 


The  dissipation  profile  obtained  by  this  method  is  designed  to  give  the 
right  shear  stress.  All  the  profile  data  mentioned  above  are  interpola¬ 
ted  with  spline  fits. 

Although  there  is  no  need  for  an  Initial  V  profile,  the  program  re¬ 
quires  initial  guesses  for  V  at  the  first  step.  Guessed  values  of  V 
should  be  provided  but  do  not  need  to  be  accurate.  Since  V  is  small, 
the  initial  guesses  were  all  taken  as  zero. 

4.3.2  Grid  Setup 

The  grid  in  the  normal  direction  at  the  initial  station  has  to  be 
provided.  Experience  shows  forty  to  fifty  grid  points  should  be  used. 
The  points  should  cover  the  region  from  y  -  0  to  approximately  y  - 
2.5  6 i ^2 •  Equal  or  variable  spacing  can  be  used.  We  found  the  best 
choice  to  be  a  compound  interest  grid: 


where  1.03  <  k  <  1.1.  This  makes  the  grid  finer  near  the  centerline 
and  coarser  near  the  edge.  In  all  the  computations  reported  herein,  45 
points  were  used  and  k  »  1.03  was  used  for  the  first  36  grids  from  the 
centerline  and  k  =  1.08  for  the  rest  of  the  domain. 

The  initial  step  size  has  to  be  specified.  A  step  size  of  about 
twice  6  j / 2  appears  to  be  reasonable.  From  the  second  step,  the  step 
size  will  be  determined  by  the  rate  the  layer  spreads  as  described 
earlier. 

4.4  Validation  of  the  Numerical  Method 


Some  validation  tests  were  performed  to  assess  the  accuracy  of  the 
present  numerical  method.  First  a  laminar  plane  jet  was  computed  with  a 
similarity  initial  condition.  The  results  were  checked  against  the  ana- 
lytlcal  similarity  solution.  With  the  45-point  compound  Interest  grid 


described  above  and  an  Initial  step  size  of  about  2  6^2*  error  is 
0.003%  in  the  centerline  streamwise  velocity  and  0. 16%  in  the  normal 
velocity  at  the  edge  after  30  initial  ^1/2*  Similar  results  were 
obtained  for  a  laminar  axisymmetric  jet. 

A  turbulent  plane  jet  calculation  was  made  to  test  the  method  when 
all  four  equations  (the  continuity,  momentum,  K  and  e  equations)  are 
coupled  and  have  to  be  solved  simultaneously.  The  results  were  compared 
with  the  similarity  solution  solved  by  Paullay  et  al.  (1985).  The 
errors  in  thickness  growth  and  centerline  velocity  decay  rates  are  both 
within  0.1%  of  the  exact  solution,  while  the  error  in  the  centerline 
value  of  K  is  about  0.9%.  When  the  number  of  grid  points  is  increased 
to  60,  almost  Identical  results  were  obtained.  This  indicates  that  grid 
independence  was  achieved  at  about  45  grid  points  across  the  layer.  A 
number  of  streamwise  step  sizes  were  tested  in  the  computation;  there  is 
no  noticeable  difference  when  the  step  size  is  26^/2  or  less. 
However,  the  difference  becomes  large  when  the  step  size  is  36j/2  or 
larger.  Therefore,  it  was  determined  that  a  step  size  of  about  1&\/2 
is  the  optimal  choice. 

4.5  Conclusions 

In  this  chapter,  an  improved  Keller  Box  method  for  solving  equa¬ 
tions  of  turbulent  free  shear  flows  has  been  described.  Several  im¬ 
provements  are  incorporated  in  the  present  method,  including  a  staggered 
grid,  an  initial  averaging  process,  and  upwind  differencing  near  the 
edge  of  the  layer.  A  grid  that  expands  with  the  layer  is  used.  A  vari¬ 
able  step  size  is  determined  by  the  growth  of  the  layer.  A  validation 
process  was  performed  to  verify  the  accuracy  of  the  method  and  the  pro¬ 
gramming. 

Initial  conditions  for  U,  K,  and  e,  an  initial  guess  for  V 
and  the  grid  at  the  initial  station  are  required,  and  their  construction 
was  described. 


Chapter  V 

ZONAL  MODELING  FOR  PLANE  FREE  SHEAR  FLOWS 


5. 1  Int  roduct Ion 

In  Chapter  III,  we  examined  experimental  data  on  the  flows  to  be 
studied  and  identified  useful  data.  In  this  chapter,  with  these  data  as 
a  guide,  we  shall  try  to  better  understand  the  physics  of  free  shear 
flows  and  to  construct  improved  models  for  these  flows  within  the  zonal 
modeling  framework.  Our  goal  is  to  build  a  consistent  zonal  model  which 
can  be  used  for  all  free  shear  flows.  However,  since  axisymmetric  flows 
differ  so  much  from  plane  flows,  we  shall  focus  our  attention  on  plane 
flows  in  this  chapter.  After  the  axisymmetric  flows  are  considered  in 
next  chapter,  we  shall  bring  these  two  cases  together. 


When  there  is  only  one  important  component  of  the  velocity-gradient 

tensor  and  the  turbulent  kinetic  energy  production  and  dissipation  rates 

are  approximately  In  balance,  the  standard  K-e  model  usually  gives 

acceptable  predictions  (Launder  et  al . ,  1972).  The  plane  pure  jet  is 

one  example.  However,  when  the  shear  is  weak,  the  standard  K-e  model 

predicts  too  slow  a  spreading  rate  of  free  shear  flows.  To  remedy  this, 

in  his  "extended"  K-e  model,  Rodi  (1972)  made  the  model  constant  Cp  a 

function  of  the  ratio  of  rates  of  turbulence  production  and  dissipation. 

This  improves  the  predictions  of  the  global  parameters  of  weak  shear 

flows  considerably.  However,  the  predicted  profiles  do  not  agree  with 

the  experiments.  Patel  and  Scheuerer  (1982),  in  their  prediction  of  a 

far  wake,  used  a  corrected  eddy  viscosity  with  an  intermittency  factor. 

This  model  improves  the  velocity  profile  but  underpredicts  the  shear 

stress  and  hence  the  spreading  rate.  Another  difficulty  with  this  model 

is  that  the  intermittency  factor  has  to  be  obtained  from  experiments  and 

is  different  from  flow  to  flow.  Others  selected  the  model  constants 

C  or  C  as  targets  for  modification  in  various  flows.  If  done  on 
e  1  e  2 

a  flow-by-flow  basis,  this  process  removes  the  possibility  of  obtaining 
a  model  valid  for  all  flows. 


It  seems,  at  the  present  time,  that  there  is  no  single  model  that 
can  predict  all  the  plane  free  shear  flows  equally  well.  This  provides 
motivation  for  the  present  study. 
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5.2  Analysis  and  Physics 


Consider  a  steady  symmetric  plane  jet  in  a  uniform  stream.  With 
the  thin  shear-layer  approximation,  the  governing  equations  for  the  flow 
are  identical  to  those  for  boundary  layers.  Beyond  the  potential  core, 
the  development  of  the  flow  is  independent  of  the  initial  conditions, 
and  the  profile  of  the  jet  becomes  self-similar,  i.e., 

U  =  U  +  U.  f (n)  (5.2.1) 

«  lm 

where 


(5.2.2) 


*1/2 


In  this  flow,  the  effect  of  the  initial  conditions  persists  only  a  rela¬ 
tively  short  distance  and  can  be  accounted  for  by  use  of  an  effective 


origin  of  the  flow  xQ.  Assuming  self-similarity,  Eqs.  (5.2.1)  and 


(5.2.2),  we  can  write  from  dimensional  analysis  (see  Bilger,  1968): 

dU, 


lm 


dx 


where 


I 


u(u  -  U  )  dy 


J  is  the  excess  momentum  flux  in  the  jet.  This  result  can  also  be  de¬ 
rived  from  the  momentum  integral  equation.  In  nondimenslonal  form,  it 
becomes 


dTITe)  '  f2(u-/uJ 


or 


lm 


(5.2.3) 


where  0  is  defined  in  Eq.  (3.2.7)  and  Is  constant  throughout  the  flow. 
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By  similar  reasoning,  we  also  obtain 


<Sl/2 

0 


(5.2.4) 


Here,  F  and  G  are  universal  functions  for  all  jet  flows.  This 
result  is  confirmed  by  the  experiment  of  Bradbury  and  Riley  (1967). 

If  the  jet  self-preserving,  the  velocity  can  be  expressed  by 
Eq.  (5.2.1)  and  shear  stress  distribution  must  have  the  form: 


-u'  v' 


u,  g(n) 

lm 


(5.2.5) 


Substitution  of  these  similarity  expressions  into  continuity  and  momen 
turn  equations  and  elimination  of  the  transverse  velocity  V  result  in: 


6l/2  dUlm 


lm 


dx 


f(e+f) 


dS 


1/2 


dx 


Bnf' 


-(- 

vu 


6l/2  dUlm 


lm 


dx 


fdq 


(5.2.6) 

where  g  is  defined  by  Eq.  (5.2.5)  and  f  by  Eq.  (5.2.1).  From  the 
momentum  integral  equation,  we  have  the  constraint: 


6 1/2  dUlm  =  d<Sl/2 

U.  dx  ^  dx 

lm 


(5.2.7) 


where  q  is  a  function  of  &'• 


Here,  1^  and  I2 
then  he  written  as 


are  defined  by  Eq 


+  IjB 


+  IjB 


(3.2.11). 


(5.2.8) 

Equation  (5.2,6)  can 


(5.2.9) 

Since  6^/2  and  8  are  functions  of  x,  the  only  conditions  under 
which  self-similarity  is  obtained,  i.e.,  f  and  g  are  independent 
of  x,  are  (see  Townsend,  1956,  and  Hill,  1965): 


d<5 


1/2 


dx 


q£3f  +  qf  -  Bnf  ’  “  (q+1 )  f 


■  f 


8  “  0  and  8  +  “  , 
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which  correspond  to  a  "pure"  jet  and  a  far  wake,  respectively.  For 
values  of  8  other  than  these  limits,  g  changes  as  8  varies.  On 
the  other  hand,  f  seems  to  be  very  insensitive  to  8,  as  shown  by 
experiments.  Moreover,  since  F  and  G  in  Eqs.  (5.2.3)  and  (5.2.4) 
are  universal  functions,  8  and  &\/2  are  related;  i.e.,  for  a  given 
8  there  is  only  one  corresponding  &\/2  an<*  v*ce  versa.  Therefore, 
the  shear  stress  level,  g,  is  only  a  function  of  6. 

As  8  increases,  the  experimental  data  show  that  g  increases  and 
the  n~value  corresponding  to  maximum  g  moves  outwards.  This  indicates 
that  the  turbulence  decays  more  slowly  than  the  centerline  velocity.  It 
also  suggests  that  diffusion  is  increased  as  8  increases;  this  is  re¬ 
quired  to  allow  the  maximum  shear  stress  location  to  move  outwards 
faster  than  the  jet  or  wake  spreads.  It  will  be  helpful  to  consider  the 
physical  reasons  why  increasing  8  has  these  effects. 

Physically,  8  represents  the  ratio  of  the  global  inertial  force 
to  the  global  driving  force  for  creating  shear  of  the  free  shear  flow. 
Thus,  when  8  increases,  the  large  eddies  become  relatively  less  ener¬ 
getic  compared  to  the  inertia  of  the  free  stream,  and  it  is  harder  for 
them  to  penetrate  into  the  free  stream.  As  a  result,  the  entrainment 
rate  Is  reduced  and  the  flow  becomes  more  parallel.  In  the  limit 
8  +  «,  the  small  eddies  are  more  responsible  for  the  turbulence  trans¬ 
port,  i.e.,  diffusion  Is  more  important.  Alternatively,  8  can  be 
regarded  as  a  ratio  of  eddy  turn-over  time  to  the  characteristic  free 
stream  time.  As  this  ratio  increases,  the  mixing  process  is  slower 
relative  to  the  free  stream  motion,  and  the  flow  spreads  less  rapidly. 
The  slow  mixing  at  high  values  of  8  also  makes  the  flow  more  intermit¬ 
tent.  It  is  well  known  that  intermittent  flow  exists  only  near  the  edge 
of  a  "pure"  jet,  whereas  it  penetrates  almost  to  the  centerline  of  a  far 
wake  (Townsend,  1956).  Interraittency  increases  the  gradient  of  the 
average  turbulence  intensity,  and  therefore  the  average  diffusion  rate 
is  increased. 

As  pointed  out  in  Chapter  III,  it  has  been  observed  by  many  re¬ 
searchers,  including  Townsend  (1956)  and  Bradbury  and  Riley  (1967),  that 
in  the  far  field,  a  co-flowing  jet  behaves  essentially  the  same  as  a 
wake.  This  means  that  the  turbulence  transport  is  independent  of  the 


sign  of  Che  mean  velocity  gradient.  Intuitively,  we  say  that  the  turbu¬ 
lence  rides  on  the  mean  flow  and  only  "sees"  a  local  velocity  difference 
without  "realizing"  its  direction.  Although  the  eddies  turn  in  opposite 
directions,  the  turbulence  transport  mechanism  appears  to  be  the  same 
for  both  cases. 

From  these  arguments,  we  conclude  that  the  turbulence  of  free  shear 
flows  is  characterized  not  by  whether  the  flow  is  a  jet  or  a  wake,  but 
rather  by  the  value  of  0.  Sabin  (1965)  also  found  a  strong  effect  of 
0  on  mixing  rate  over  a  wide  range  of  values  of  0  for  the  mixing 
layer.  Using  similarity  arguments,  he  provided  a  general  relation  be¬ 
tween  spreading  rate  and  0.  The  same  relation  was  found  by  Abramovich 
(1963),  apparently  independently.  This  Sabin-Abramovich  relation  was 
found  to  be  one  of  the  few  universally  valid  correlations  in  both  the 
1972  NASA-Langley  Conference  on  Evaluation  of  Computation  in  free  shear 
layers  and  the  1980-81  AFOSR-Stanf ord  Conference  evaluating  complex  flow 
fields.  Thus  the  importance  of  0  (or  some  wholly  equivalent  parameter 
of  different  mathematic  form)  is  fully  validated  and  very  widely  accep¬ 
ted  for  the  mixing  layer  case. 

All  these  considerations  suggest  that  a  turbulence  model  designed 
to  simulate  these  flows  in  a  consistent  way  should  depend  on  0  and  not 
the  mean  flow  type. 

In  summary,  we  note  the  following: 

1.  The  turbulence  model  should  reflect  the  relative  importance  of 
the  diffusion  process  as  indicated  by  the  value  of  0. 

2.  It  appears  that  we  can  consider  regions  with  the  limits  of  0 

0  and  0  ♦  =°  as  two  zones  of  free  shear  flows  and  intermedi¬ 
ate  values  of  0  as  a  readjustment  zone. 

3.  A  co-flowing  jet  or  a  wake  starts  with  a  small  value  of  0  and 
develops  a  larger  0  downstream.  In  terms  of  our  model,  it 
starts  in  the  first  zone  and  readjusts  toward  the  second.  A 
"pure"  jet,  on  the  other  hand,  stays  in  the  first  zone  perma¬ 
nently.  A  given  mixing  layer  has  a  constant  value  of  0,  but 
mixing  layers  as  a  class  have  values  of  0  running  from  1  to 

To  apply  the  zonal  modeling  approach,  we  shall  model  the 
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two  limiting  cases  first  and  then  patch  them  together  by  a 
blending  technique  which  describes  the  readjustment. 


5.3  Construction  of  Zonal  Models 

As  noted  in  Chapter  I,  the  standard  K-e  model  will  be  used  as  a 
base  for  the  development  of  zonal  models.  Therefore,  in  constructing 
zonal  models  we  shall  first  test  the  standard  K-e  model  in  each  case. 
If  it  works  for  a  given  case,  no  modifications  will  be  made,  since  it  is 
then  a  satisfactory  zonal  model. 


5.3.1  Zonal  Model  for  the  Limiting  Case  g  •»  0 

A  "pure"  jet  represents  the  paradigmatic  flow  for  this  zone.  After 
the  disappearance  of  the  potential  core,  a  pure  jet  quickly  becomes 
self-preserving.  Therefore,  we  can  use  pure  jet  data  at  any  downstream 
location  to  tune  the  zonal  model  for  this  zone. 

The  standard  K-e  model  was  tested  for  this  flow  and  gives  sur¬ 
prisingly  good  results.  The  spreading  rate  and  centerline  velocity 
decay  rate  are 


dx 


0. 108 


dx 


0.16 


which  are  in  excellent  agreement  with  Eqs.  (3.2.4)  and  (3.2.5).  The 
mean  velocity,  turbulent  kinetic  energy,  and  shear  stress  profiles  also 
agree  quite  well  with  experimental  data.  Figures  5.1  to  5.3  illustrate 
the  agreement. 

The  standard  K-e  model  is  known  to  work  well  for  flows  with  just 
one  significant  velocity  gradient  in  which  the  turbulent  kinetic  energy 
production  rate  approximately  equals  the  dissipation  rate.  The  plane 
pure  jet  satisfies  these  two  conditions.  Another  part  of  the  reason 
that  the  standard  K-e  model  predicts  pure  jet  flows  so  well  is  that 
some  of  the  model  constants  were  based  on  pure  jet  data.  There  is  no 
need  to  modify  the  standard  model  for  this  zone. 


5.3.2  Zonal  Model  for  the  Limiting  Case  B  *  00 

A  far-field  wake  or  co-flowing  jet  both  are  instances  of  this  zone. 
However,  Ujm  is  finite  at  the  last  measuring  station  in  all  experi¬ 
ments.  Fortunately,  when  Ulm  is  less  than  10?  of  U,,,  the  flow 
behaves  self-similarly  and  the  asymptotic  far-field  state  is  closely 
approximated.  Therefore,  measurements  which  satisfy  this  criterion  can 
be  used.  In  the  co-flowing  jet  of  Bradbury  and  Riley  (  1967),  Ujm  is 
never  less  than  20%  of  U^;  this  is  not  small  enough  to  qualify.  Pot's 
(  1979)  wake  data,  on  the  other  hand,  reach  0^/U^  “  5%.  Therefore,  the 
measurements  at  the  last  few  stations  in  Pot's  experiment  will  be  used 
to  tune  the  model  for  this  zone.  In  this  range,  the  global  parameters 
of  Pot's  data  behave  as  described  by  Eqs.  (3.2.15)  and  (3.2.16). 


The  standard  K-e  model  was  first  tested.  The  results  are: 


d(s1/2/e)2 

d(x/9) 


0.065 


d(x7e) 


0.248 


These  quantities  are  significantly  lower  than  those  of  Eqs.  (3.2.15)  and 

(3.2.16) ,  respectively.  The  spreading  parameter,  as  defined  in  Eq. 

(3.2.17) , 


sp  -  0.0635 

is  more  than  35%  lower  than  Pot's  and  Townsend's  values.  In  addition, 
the  predicted  mean  velocity,  turbulent  kinetic  energy  and  shear  stress 
profiles  do  not  agree  with  the  data.  Figures  5.4  to  5.6  show  the  com¬ 
parisons.  Not  only  are  the  magnitudes  in  serious  error,  the  shapes  are 
also  not  correct,  especially  near  the  edge  of  the  layer.  These  various 
profiLes  all  go  to  zero  too  fast. 

From  the  discrepancies  between  the  K-e  predictions  and  the  exper¬ 
iments,  it  appears  that  the  diffusion  predicted  by  the  standard  K-e 
model  is  not  sufficiently  strong  for  this  zone.  To  remedy  this,  we 
decided  to  decrease  two  model  constants:  C  and  Ou .  These  changes 

t.  J 

increase  the  shear  stress  as  well  as  diffusion.  Although  this  choice  of 


model  constants  may  not  be  unique,  It  seems  the  simplest  way  to  achieve 
the  necessary  effect.  This  satisfies  our  guideline  that  we  want  to  make 
the  changes  as  simple  as  possible  when  modification  is  necessary. 


After  a  series  of  tests  against  the  last  few  stations  of  Pot's  far 
wake  data,  we  set  the  two  model  constants  at: 


and 


1.04 


(5.3.  1) 


0.65 


(5.3.2) 


Using  these  values,  all  predicted  profiles  agree  reasonably  well  with 
the  experimental  data  as  shown  in  Figs.  5.4  to  5.6.  Moreover,  the 
global  parameters  match  those  described  by  Eqs.  (3.2.15)  and  (3.2.16) 
within  3%.  They  are: 

d(<Sl/2/e)2 


d(x7e)" 


0. 103 


d(U./Ulm) 

dYx/e ) 


-  0.437 


Thus,  we  have  constructed  a  model  for  the  zone  8  ♦  »  by  changing 
two  model  constants. 

5.3.3  Model  for  the  Readjustment  Region 

We  have  now  zonal  models  for  the  limiting  cases  6*0  and  8  *  •• 
Next,  we  need  to  combine  these  to  form  a  complete  model.  Any  flow  with 
8  between  zero  and  infinity  can  be  regarded  as  being  in  a  transitional 
stage  between  the  two  zones.  This  transition  is  called  a  "readjustment" 
in  zonal  modeling.  To  reflect  the  continuous  change  of  physics  in  read¬ 
justment,  we  shall  patch  the  two  zonal  models  such  that  the  model  con¬ 
stants  vary  smoothly  between  the  limiting  values.  The  experience  with 
homogeneous  flows  suggests  that  the  two  constants  can  be  expressed  as: 


1.44  -  0.4  X 


1 


-  1.0  -  0.35  Xj 
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(5.3.3) 

(5.3.4) 


where  0  _<_  A  j  1  is  the  readjustment  parameter.  Aj  is  a  function 
of  3  and,  again  in  line  with  earlier  experience,  we  chose  it  to  be  the 
solution  of  a  first-order  differential  equation: 


dA, 

dT  ‘  c<‘  -  V 


(5.3.5) 


where  C  is  a  constant  to  be  determined.  The  solution  to  Eq .  (5.3.5) 


A 1  *  1 


(5.3.6) 


Note  that  the  model  constants,  Eqs.  (5.3.3)  and  (5.3.4),  reduce  to  the 
appropriate  values  for  each  of  the  two  zones  as  3  approaches  the  lim¬ 
iting  values. 

The  exponential  form  of  Eq.  (5.3.6)  makes  the  constants  leave  the 
first  zone  (8  ♦  0)  values  rather  rapidly  and  reach  the  second  zone 
values  asymptotically.  This  seems  to  agree  with  the  data.  It  is  veri¬ 
fied  by  the  rapid  growth  of  the  maximum  turbulent  kinetic  energy  and 
shear  stress  in  the  near  field,  as  shown  in  Figs.  5.9  and  5.10. 

The  data  at  intermediate  8  values  (8  <  10  or  U,  /U  >0.1) 

in  00 

in  Pot's  wake  experiment  are  excellent  for  tuning  the  readjustment  con¬ 
stant  0.  With  the  use  of  these  data,  it  was  found  that 

C  -  0.18  (5.3.7) 


This  value  of  C  gives  good  agreement  with  the  data  and  completes  the 
model  for  planar  free  shear  flows. 

Predictions  of  the  zonal  model  are  compared  with  Pot's  wak«  data, 
for  the  spreading  rate,  centerline-velocity-deficit  decay  rate,  the 
maximur  kinetic  energy,  and  the  maximum  shear  stress  in  the  layer, 
respectively,  in  Figs.  5.7  and  5.10.  The  agreement  is  excellent.  Also 
presented  in  these  figures  are  the  results  predicted  by  the  standard 
K-e  model.  The  disagreement  between  the  data  and  the  standard  K-e 
model  is  significant,  especially  in  the  far  field. 


5 . 4  Tests  of  the  Zonal  Model 

We  have  constructed  a  complete  model  by  patching  two  zonal  models 
for  the  extreme  cases.  This  model  was  tuned  using  only  Pot's  wake  data. 
This  model  gives  very  good  agreement  with  that  data  set  as  seen  in  the 
preceding  section.  However,  the  purpose  of  modeling  is  to  enable  us  not 
only  to  utilize  the  model  for  flows  to  which  it  is  tuned,  but  also  to 
other  flows.  We  use  co-flowing  jets  and  mixing  layers  ror  this  check  in 
this  sec t i on . 

5.4.1  Co-flowing  Jets 

It  was  emphasized  in  Chapter  III  that  behavior  of  co-flowing  jets 
may  vary  from  "pure"  jet  flow  to  self-preserving  wake  flow.  This  flow 
provides  a  good  test  for  the  present  model.  The  data  set  recorded  by 
Bradbury  and  Riley  (1967)  was  reviewed  in  Chapter  III  and  seen  to  be 
consistent  and  of  high  quality.  Also,  the  data  were  nondimensional ized 
to  collapse  all  co-flowing  Jet  data  with  different  ratios  of  nozzle 
exhaust  to  free  stream  velocity.  We  shall  use  this  data  set  in  the 
test . 

An  Initial  condition  was  generated  using  the  data  of  Bradbury's 
earlier  experiment  (1965).  Calculation  was  made  to  a  downstream 
location  of  about  100  0.  The  results  for  the  Jet  spreading  rate, 
centerline  velocity  decay,  and  the  eddy  Reynolds  number  defined  by  Eq . 
(5.2.14)  are  shown  in  Figs.  5.11  to  5.13,  respectively.  The  predictions 
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pared.  The  standard  K-e  model,  on  the  other  hand,  predicts  well  only 


In  the  Initial  part  of  the  flow  and  departs  from  the  data  further  down¬ 
stream.  This  is  expected  because,  as  demonstrated  earlier,  the  stan¬ 
dard  K-e  model  works  well  for  flows  with  small  8  (the  first  zone) 
but  poorly  for  flows  with  large  8.  The  serious  disagreement  between 
the  data  and  the  standard  K-t  model  prediction  for  the  eddv  Reynolds 


number  seen  in  Fig.  5.13  reveals  that  the  standard  model  is  unaMe  to 
predict  the  shear  stress  level  accurately.  The  low  shear  stress  pro¬ 
duces  the  slower  spreading  and  velocity  decay. 


Another  challenging  test  is  the  mixing  layer.  These  flows  do  not 
have  similar  velocity  profiles.  However,  the  ratio  of  inertial  force  to 
driving  force  can  be  defined  if  we  consider  to  be  the  average  of 
the  two  free  stream  velocities: 


{  (U.  +  ) 


(5.4.1) 


and  U,  to  be  the  difference  between  either  free  stream  velocity  and 
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Thus,  we  take  8  to  be  U  /U.  -  (U  +U  )/(U„  -U« 111  these  expres_ 
slons,  and  are  the  free  stream  velocities  on  the  high  and 
low  velocity  sides,  respectively.  One  side  has  a  velocity  excess;  the 
other  has  a  velocity  defect.  Both  sides  share  the  common  and  have 
the  same  value  of  Uj^,  and  hence  identical  8's. 

Another  way  of  viewing  a  mixing  layer  is  as  follows.  There  exists 
a  dividing  streamline  which,  if  not  exactly  coincident  with,  is  very 
close  to  the  line  with  velocity  equal  to  l'^  .  The  high  velocity  side 
can  be  regarded  as  jet-like  and  the  low  velocity  side  as  wake-like. 

Using  the  velocities  as  defined  in  Eqs.  (5.4.1)  and  (5.4.2),  the 
ratio  of  Inertial  force  to  driving  force,  8,  for  the  mixing  layer  Is 
Identical  to  the  parameter  used  by  Sabin  (1965)  to  correlate  the  spread¬ 
ing  parameter  tor  Mows  with  various  velocity  ratios. 

In  the  calculation  of  mixing  layers,  we  have  to  note  that  the 
results  are  verv  sensitive  to  the  boundary  conditions  on  the  normal  vel¬ 
ocity,  V,  at  both  the  top  and  bottom  boundaries,  as  explained  in  Chap¬ 
ter  III.  In  ’rder  to  aciurately  compute  these  flows,  therefore,  both 
boundary  conditions  should  be  specified.  However,  in  the  thin  shear- 
!aye>-  ipnrox  i  ma  t  1  on ,  only  one  V  boundary  condition  can  be  given.  There 
is  no  wav  to  Improve  the  approximation  other  than  to  select  a  V  bound¬ 
ary  condition  that  is  best  for  a  particular  problem.  There  a-e  several 
wavs  to  specify  this  boundary  condition,  Including 
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In  the  above  expressions,  V  and  V_  are  the  normal  velocities  on 
r  °*1  "2 

the  top  and  bottom  boundaries,  respectively.  The  last  is  derived  from 
the  y-momentum  conservation  equation  and  is  the  best  condition  for  a 
"free"  mixing  layer,  l.e.,  one  that  is  not  affected  by  wind-tunnel 
walls. 


The  V  boundary  condition  determines  the  angle  at  which  the  mixing 
layer  propagates  and  can  be  chosen  to  match  the  calculated  flow  angle  to 
that  observed  experimentally;  this  differs  from  experiment  to  experi¬ 
ment.  Fortunately,  the  spreading  rate  is  only  slightly  affected  by 
the  V  boundary  conditions.  In  this  work,  only  the  spreading  rate  is 
of  major  concern;  the  flow  angle  will  not  be  examined. 

Mixing  layers  with  various  velocity  ratios  are  computed  with  both 
the  present  model  and  the  standard  K-e  model.  Calculations  were 
carried  out  until  a  self-similar  solution  was  achieved.  Figure  5.14 
shows  the  comparison  between  the  predicted  spreading  rates  and  those  of 
various  experiments  for  the  range  of  1  <  8  <  9.  The  present  model 
performs  better  than  the  standard  K-e  model.  Figures  5.15(a),  (b)  and 
(c)  compare  the  shear  stress  profiles  of  the  model  predictions  and  the 
data  for  three  different  velocity  ratios  (R  ”  0,  0.3,  and  0.61), 
respect  1 vely . 


5 . 5  Galilean  Invariance  —  A  Frame-Invariant  Form  of  the  Model 

The  description  of  a  phenomenon  must  be  the  same  in  all  inertial 
frames  of  reference.  That  is,  the  equations  governing  the  physics 
should  be  Identical  whether  the  reference  frame  Is  stationary  or  is 
moving  at  a  constant  speed.  To  be  more  specific,  the  equations  must  be 
form-invariant  under  the  Galilean  tranf ormat ion 


£  - 


(5.5.  1  ) 


where  £ '  and  X  are  the  coordinates  in  the  new  and  old  reference 
frames,  respectively.  The  new  reference  frame  has  a  constant  relative 
velocity  $  with  respect  to  the  old  frame.  The  Navier-Stok.es  equa¬ 
tions,  for  example,  are  Galilean  invariant  in  this  sense. 

It  is  essential  that  the  turbulence  model  be  Galilean-invariant. 
Otherwise,  the  model  will  give  different  results  in  various  inertial 
reference  frames.  The  model  developed  for  plane  free  shear  flows  in 
this  chapter  Is  not  Ga 1 i lean- i nva r i an t  because  the  parameter  6  on 
which  the  model  constants  depend  contains  U  ,  which  changes  with  the 


J 


frame  of  reference.  This  would  restrict  use  of  the  model  to  reference 
frames  in  which  the  flow  appears  steady.  For  example,  for  jet  flows, 
the  reference  frame  has  to  be  stationary  with  respect  to  the  nozzle  that 
produces  the  jet  and,  for  wake  flows,  stationary  with  respect  to  the 
hodv  that  generates  the  wake.  For  all  other  reference  frames,  the  pre¬ 
sent  model  would  be  invalid. 

A  look  at  the  analysis  indicates  that  the  lack  of  Galilean  Invar¬ 
iance  of  the  node!  arises  from  our  assumption  of  a  steady  flow.  When 
•he  reference  frame  is  changed,  the  flow  appears  unsteady  and  the 
analysis  of  Sec.  ^ .  2  is  no  longer  adequate.  Therefore,  in  order  to 
obtain  Galilean  invariance,  we  mils'  consider  unsteady  flows. 

The  momentum  integral  equation  for  an  unsteady  flow  is' 
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Combining  Eqs.  (5.5.2)  and  (5.5.3)  gives 


4  1/2  3L’h 


f(B  +  f) 


Sl/2  3U1, 
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which  is  identical  to  Eq .  (5.2.3).  Moreover,  the  parameter  B  defined 
by  (5.5.3)  reduces  to  the  previous  definition,  Eq .  (5.2.1),  for  a  steady 
flow.  The  remainder  of  the  analysis  is  the  same  as  that  in  Sec.  5.2  and 
will  not  be  repeated  here. 

Therefore,  in  applying  the  present  model  to  a  general  free  shear 
flow,  the  definition  of  the  parameter  8  given  by  Eq .  (5.5.3)  has  to 

replace  that  of  F.q .  (5.2.1)  This  aak.es  the  model  Ca  1 1 1  ean-i  nvar  i  ant . 

No  other  changes  are  needed.  Since  the  original  equations  with  8 
given  bv  Eq.  (5.2.1)  are  simpler,  the  more  complex  analysis  of  this 
section  need  be  used  only  when  the  flow  is  unsteady  in  the  reference 
f  rame  chosen . 

5 . b  Cone  1  us  1 ons 

In  this  chapter,  plane  free  shear  flows  were  studied.  A  zonal 
model  was  developed.  Two  zones  were  identified;  they  can  be  called  the 
■.ear  and  far  fields,  although  that  nomenclature  is  not  accurate  for  all 
Mows.  This  problem  is  a  semantic  one  arising  from  the  definition  of  a 
zone,  and  we  shall  suggest  a  clearer  semantics  in  the  final  chapter. 
The  difference  in  the  turbulence  transport  mechanism  in  the  two  zones 
was  discussed.  Zonal  models  have  been  constructed  for  these  two  zones. 
A  patching  technique  was  used  to  represent  the  readjustment  and  to  link 
the  two  zonal  models.  The  new  model  predicts  the  change  in  physics  by 
hanging  two  constants.  The  requirement  of  Calilean  invariance  can  be 
satisfied  by  use  of  a  more  general  definition  for  8.  Test  results  show 
that  the  model  developed  in  this  chapter  performs  much  better  than  the 
standard  K-t  model,  but  is  no  more  difficult  to  apply  or  expensive  to 
run  for  plane  free  shear  flows. 
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Chapter  VI 


ZONAL  MODELING  FOR  AXISYMMETRIC  FREE  SHEAR  FLOWS 


b  .  1  Int roduct  ton 

In  Chapter  V,  we  built  a  new  (zonal)  mode L  that  works  well  for  all 
usable  cases  of  plane  free  shear  flows.  This  new  model  Improves  the 
standard  K-e  model  significantly  by  allowing  two  of  the  constants  to 
become  functions  of  the  parameter  6.  Axisymmet ric  flows  were  not  con¬ 
sidered  in  Chapter  V,  because  there  is  a  fundamental  difference  in  the 
turbulence  structure  due  to  the  lateral  vortex  stretching  in  axisymmet- 
ric  Mows  as  mentioned  in  Chapter  III  (Fig.  3.12).  When  vortex  stretch¬ 
ing  is  significant,  it  decreases  the  ratio  of  'u'v'i/K  which  deter¬ 
mines  a  central  constant  in  the  K-e  model.  In  this  chapter,  we  shall 
examine  the  vortex  stretching  effect  in  detail  and  incorporate  its 
effe'ts  Into  the  model  for  axisymmet ric  flows  through  a  suitable  non- 
dimensional  p  l  r  a  me  t  e  r  . 

The  standard  K-e  model  was  tuned  to  plane  shear  flows.  As  a 
result,  boundary  layers  and  plane  pure  Jets  are  well  represented  by  this 
model.  However,  to  predict  ax i sv mme t r l c  flows,  some  modifications  are 
required.  Yost  of  the  modifications  shown  in  the  literature  involve 

hanging  either  C  or  C  ror  example,  Launder  et  al .  (1972)  added 
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No  convincing  physical  explanation  is  provided  to  justify  any  of 
these  nod l f icat ions .  Moreover,  when  these  corrections  are  used,  the 
generality  of  the  model  ability  to  predict  plane  and  axi symmet r i c  flows 
is  destroyed.  It  is  not  known  whether  the  correction  should  be  used  for 
flows  Intermediate  between  plane  and  axisymmet r ic .  Moreover,  these 
modifications  were  primarily  tuned  to  pure  jets;  they  do  not  predict 


wakes  well.  In  his  extended  K-e  model,  in  which  is  a  function 


of  the  ratio  of  the  turbulence  production  and  dissipation  rates,  Rodi 
(1972)  included  a  modification  for  axisymmetric  flows  similar  to  that  of 
l.aunder  at  el.  (1972)  described  above.  This  gives  predictions  for  axi- 
svmmetri-  flows  better  than  the  ones  mentioned  above.  However,  the  pre¬ 
dicted  profiles  do  not  agree  well  with  data. 


In  this  chapter,  we  will  use  the  zonal  approach  to  develop  a  model 
for  axisymmetric  flows.  We  want  to  maintain  the  ability  to  handle  both 
'lane  and  axisymmetric  cases.  Therefore,  we  shall  often  refer  to  the 
plane  flow  model  developed  in  Chapter  V.  Our  objective  is  to  create  a 
general  mode  1  which  treats  plane  flows  as  one  special  case  and  axisym¬ 
metric  flows  as  another. 


h  . Ann  1  v sis  and  Physics 

An  analysis  similar  to  the  one  for  plane  flows  (Sec.  5.2)  can  be 
carried  out  for  axisymme t ri c  flows.  To  briefly  illustrate  several  key 
points,  consider  a  steady  axisymmetric  jet  in  a  uniform  stream  uith 
arbitrary  free  stream  velocity.  i’sing  dimensional  analysis,  we  find 
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Here,  0  is  the  momentum  thickness  for  axisymmetric  flows  defined  by 
Eq .  (3.3.4).  F  and  C  are  universal  functions  for  all  axisymmetric  jets. 
Indeed,  Rodi  (1975),  using  Eqs.  (6.2.1)  and  (6.2.2),  collapsed  the  jet 
width  and  centerline  velocity  data  of  three  experiments  (Fig.  3.13). 

We  can  apply  the  similarity  expressions  used  in  Chapter  V  for 
streamwise  velocity  and  turbulent  shear  stress.  Substituting  them  into 
the  momentum  equation,  we  have 
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The  momentum  integral  equation,  after  the  similarity  transformation, 
reads 
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and  I{  and  are  defined  by  Eq.  (3.3.6).  Combining  Eqs.  (6.2.3) 

and  (6.2.4)  results  In 
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Since  B  and  5^2  are  related  by  Eqs.  (6.2.1)  and  (6.2.2),  the  shear 
stress  can  be  considered  a  function  of  B  alone. 
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Strictly  speaking,  the  parameter  B  =  Uoo/Uim  can  be  used  or,ly  in 

steady  flows.  In  general,  g  =  (D5 , /Dt ) /(U ,  36,._/3x)  should  be 

1/2  lm  1/2 

used  to  allow  Galilean  invariance.  This  point  was  discussed  in  the  last 
chapter. 


The  axisymmetric  wake  can  be  analyzed  in  a  similar  manner.  Let  us 
assume  that  the  turbulence  transport  mechanism  is  independent  of  the 
sign  of  mean  velocity  gradient  and  only  a  function  of  6*  As  in  plane 
flows,  0  has  two  extreme  values:  zero  and  infinity.  These  correspond 
to  a  pure  jet  and  a  far  wake,  respectively.  A  co-flowing  jet  or  wake 
will  develop  between  these  two  limits.  They  may  start  from  a  state 
close  to  a  pure  jet  (small  0)  and  eventually  arrive  at  the  asymptotic 
far  wake  state  (large  0).  Therefore,  the  zonal  idea  used  in  plane 
flows  can  as  well  be  applied  here.  We  shall  first  consider  these  two 
extremes  (or  zones). 

The  first  zone  (0  *  0)  represents  a  pure  jet.  This  flow  spreads 

linearly.  Its  spreading  rate  is  the  greatest  among  axisymmetric  flows. 
Due  to  the  rapid  spreading  and  small  radius  of  curvature  of  the  flow, 
there  exists  a  significant  strain  in  the  circumferential  direction. 
This  strain  tends  to  stretch  the  vortex  rings  which  form  in  the  shear 
layer.  For  this  effect  to  be  important,  two  criteria  must  be  met. 
Firstly,  the  jet  has  to  have  a  radius  that  is  not  large  in  comparison 
with  the  shear  layer  width.  Secondly,  the  flow  has  to  spread  rapidly. 
Otherwise,  the  straining  will  be  weak. 

The  stretching  of  vortex  rings  reduces  the  ratio  of  |u'v'|/K. 
This  conclusion  is  reached  on  the  basis  of  experimental  data.  No  physi¬ 
cal  or  mathematical  derivation  is  available  at  the  present  time.  For  a 
given  turbulence  intensity,  the  vortex  stretching  causes  the  shear 
stress  to  decrease.  For  axisymmetric  pure  jets,  in  which  the  effect  of 
vortex  stretching  is  the  strongest,  the  ratio  of  |u'v'|/K  in  the 
equilibrium  region  is  reduced  to  about  0.23  (vs.  0.3  for  the  plane 
flows,  cf.  Fig.  3.12). 

The  second  zone  (0  +  °°)  is  an  asymptotic  far  wake.  In  this  zone 
the  spreading  rate  is  very  slow.  The  vortex  stretching  effect  is  there¬ 
fore  small.  Indeed,  in  the  review  of  Antonia  and  Bilger's  (1973)  axi- 
symnetric  co-f lowing  jet  experiment  in  Chapter  HI,  we  found  that 
|u'v'|/K  is  close  to  0.3  in  the  equilibrium  region  far  downstream.  The 
physics  of  this  far-field  zone  is  nearly  identical  to  that  of  the  plane 
asymptotic  far  wake. 

The  turbulence  diffusion  enhancement  as  0  increases  seems  to  be 


the  same  for  axisymmetric  flows  as  for  plane  flows. 


6.3  Construction  of  Zonal  Models 

Our  approach  will  again  be  to  model  the  extremes  first,  then  the 
readjustment . 

6.3.1  Zonal  Model  for  the  Limiting  Case  B  *  0 

As  before,  the  standard  K-e  model  is  first  tested  to  determine 
what  modification,  if  necessary,  is  needed.  The  results  for  the  global 
parameters  are: 
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Both  the  spreading  and  velocity  decay  rates  are  over-predicted  by  more 
than  35%  by  the  standard  K-e  model.  Mean  velocity,  turbulent  kinetic 
energy  and  shear  stress  profiles  are  shown,  respectively,  in  Figs.  6.1 
to  6.3  and  compared  with  experimental  data.  While  the  mean  velocity  and 
kinetic  energy  profiles  show  reasonable  agreement  with  the  data,  the 
shear  stress  is  significantly  over-predicted,  i.e.,  the  ratio  |u'v'|/K 
is  too  high  for  this  flow.  The  reason  is  that  C^  =  0.09  in  the  stand¬ 
ard  K-e  model  |u'v'|/K  at  appproximately  0.3  in  the  region  where 

production  is  in  balance  with  dissipation.  However,  as  emphasized 
above,  |u'v'|/K  is  reduced  by  the  effect  of  vortex  stretching  and  is 
0.23  rather  than  0.3.  To  reflect  this  value  of  |u'v'|/K  in  the  model, 
it  appears  that  the  value  of  Cy  should  be  decreased. 

From  the  shear  stress  and  eddy  viscosity  formulation,  as  described 


in  Chapter  II,  we  have 


=  C  - 
u  e 


(6.3.1) 


where  P  =  -  u'v'  3U/3y  is  the  production  rate  of  the  turbulence. 
Knowing  that  |u'v'|/K  ■  0.23  in  the  region  where  production  Is  approx¬ 
imately  equal  to  dissipation  for  axisymmetric  pure  jet,  we  set 
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which  are  in  excellent  agreement  with  Eqs.  (3.3.1)  and  (3.3.2)  that  cor¬ 
relate  the  data.  The  various  profiles  generated  by  using  ■  0.05  are 
also  presented  in  Figs.  6.1  to  6.3.  All  three  profiles  match  the  data 
very  well. 

The  zonal  model  for  this  zone  is,  therefore,  the  standard  K-e 
model  with  a  change  in  model  constant  to  account  for  the  vortex 

stretching.  No  other  modification  is  needed. 

6.3.2  Zonal  Model  for  the  Limiting  Case  8  «> 

A  far-field  wake  or  co-flowing  jet  represents  this  zone.  As  noted 

above,  due  to  the  slow  spreading,  the  effect  of  vortex  stretching  is 

i 

negligible  in  this  zone,  and  the  physics  of  the  flow  is  nearly  the  same 
as  that  of  plane  cases.  Therefore,  the  zonal  model  of  the  plane  asymp¬ 
totic  far  wake  was  adopted  for  the  axisymmetric  counterpart.  This  allows 
maintenance  of  the  generality  of  the  new  model.  No  tuning  of  model  con¬ 
stants  is  needed  or  performed  in  this  zone. 


6.3.3  Model  for  the  Readjustment  Region 

As  stated  earlier,  our  intent  is  to  produce  a  general  model  for 
free  shear  flows.  Therefore,  we  want  to  keep  the  model  for  axisymmetric 
flows  as  consistent  as  possible  with  that  for  plane  flows.  Furthermore, 
the  changes  with  6  appear  to  be  the  same  for  both  plane  and  axisym¬ 
metric  flows.  For  these  reasons,  we  shall  keep  the  variation  of  the 
model  constants  C  and  Ot,  between  the  two  zones  the  same  as  in 

C  J  1C 

plane  flows,  i.e.,  we  use  Eqs.  (5.3.3)  to  (5.3.7).  This  leaves  us  only 


one  more  linking  to  be  made  (for  the  model  constant  C^)  before  the 
zonal  model  for  axisymmetric  flows  is  completed. 

In  going  from  the  first  zone  to  the  second,  the  stretching  dimin¬ 
ishes.  The  model  constant  Cy  should  increase  from  0.05  to  0.09.  This 
transitional  stage  can  be  represented  by  sliding  between  the  two 
extreme  values.  Before  constructing  the  blending  function  for  C^>  let 
us  examine  what  controls  the  stretching  of  vortices. 


As  discussed  in  Sec.  6.2,  the  flow  has  to  be  small  in  the  normal 
direction  and  has  to  spread  rapidly  in  order  for  the  stretching  effect 
to  be  significant.  This  means  that  the  vortex  stretching  should  be  a 
function  of  the  spreading  rate,  arid  the  radius  of  curvature 
of  the  vortex  ring.  From  these  and  dimensional  arguments,  an  appropri¬ 
ate  parameter  appears  to  be: 


S 


Rl/2  1 
6l/2  d6l/2/dx 


(6.3.3) 


where  ri/2  is  radius  of  curvature  of  the  vortex  line  passing 
through  the  point  of  6  j  y  2  *  Tlle  smaller  this  parameter,  the  more 
stretching  there  is.  For  plane  flows,  S  is  infinite  and  there  is  no 
stretching.  For  axisymmetric  flows,  ri/2  =  ^1/2  anc*  ^  *s  inversely 
proportional  to  the  spreading  rate.  For  general  three-dimensional 
flows,  for  example  a  rectangular  jet,  S  is  smaller  at  the  corners  due 
to  the  small  radius  of  curvature.  This  gives  more  stretching  and  hence 
smaller  shear  stress  near  the  corners.  Consequently,  the  spreading  near 
the  corners  slows  down  and  the  flow  will  eventually  become  axisymmetric. 
This  is  at  least  qualitatively  correct. 


The  smallest  observed  value  of  S  occurs  in  the  axisymmetric  pure 
jet  and  is  11.3.  Therefore,  we  propose  that  the  sliding  function  for 
model  constant  be: 


C  =  0.05  +  0.04  X0  ,  S  >  11.3  (6.3.4) 
U  2  — 

where  0  <  X£  *  Is  the  readjustment  parameter.  Further  study  may  be 
needed  if  a  value  of  S  below  11.3  occurs  in  any  flow.  X2  is  governed 
by  the  first-order  differential  equation: 


air  ■  c80  -  x2>  <6-3-5) 


where  Cg  is  a  constant  that  needs  to  be  tuned.  The  solution  to  Eq. 
(6.3.5)  is: 

-C  (S-11.3) 

X2  =  1  -  e  S  (6.3.6) 


The  last  task  is  to  determine  the  constant  Cg .  The  tuning  is  done 
by  computing  a  co-flowing  jet  and  comparing  the  results  with  the  data  in 
Figs.  3.13  and  3.14.  An  initial  condition  was  created  by  using  Antonia 
and  Bilger's  (1973)  experimental  data  for  a  co-flowing  jet  with  the  jet 
nozzle  velocity  to  free  stream  velocity  equal  to  4.5.  It  was  determined 
that 


Cs  -  0.1  (6.3.7) 

gives  the  best  results.  The  agreement  for  the  jet  spreading  rate  and 
centerline  velocity  decay  rate  can  be  seen  in  Figs.  6.4  and  6.5.  The 
K-e  model  prediction  is  included  in  these  figures  for  comparison. 

6.4  Tests  of  the  Zonal  Model 

We  have  constructed  a  new  model  for  axisymmetric  free  shear  flows 
with  Eqs.  (5.3.3),  (5.3.4)  and  (6.3.4)  and  are  ready  to  test  the  model. 

6.4.1  Co-Flowing  Jets 

It  was  decided  to  simulate  another  axisymmetric  co-flowing  jet. 
The  initial  condition  is  based  on  the  data  of  another  co-flowing  jet  of 
Antonia  and  Bilger  (1973).  This  jet  has  a  ratio  of  nozzle  velocity  to 
free  stream  velocity  equal  to  3.  Computation  was  carried  out  to  about 
150  0  downstream.  The  spreading  rate  and  centerline  velocity  decay  are 
presented  in  Figs.  6.6  and  6.7,  respectively.  Good  agreement  with  ex¬ 
perimental  results  is  obtained  for  both  quantities.  Predictions  of  the 
standard  K-e  model  for  the  same  flow  are  also  shown  for  comparison. 
Although  giving  reasonable  spreading  rate,  the  standard  K-e  model 
fails  to  predict  the  velocity  decay  rate  accurately. 
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6.4.2  Wakes 


Axisymmet ri c  wakes  provide  a  challenge  to  the  present  model  because 
these  flows  were  not  used  In  the  determination  of  the  model  constants. 
Chevray's  (1968)  wake  experiment  Is  used  for  this  test.  The  computation 
uses  the  data  at  x/D  =  6  as  an  Initial  condition.  Comparisons  of  the 
prediction  by  the  present  model  and  experimental  data  for  the  spreading 
rate  and  decay  rate  of  centerline  velocity  deficit  are  shown  in  Figs. 
6.8  and  6.9.  Excellent  agreement  is  achieved.  The  standard  K-e  model 
underpred lets  both  quantities;  these  results  are  also  shown  in  the  fig¬ 
ures.  The  predicted  turbulent  kinetic  energy  and  shear  stress  profiles 
at  the  last  measuring  station  (x/D  =  18)  are  compared  with  the  exper¬ 
imental  data  in  Figs.  6.10  and  6.11.  Good  agreement  is  observed  with 
the  present  model.  The  standard  K-e  model,  on  the  other  hand,  pre¬ 
dicts  significantly  lower  values  for  both  quantities;  the  low  shear 
stress  is  responsible  for  the  slow  spreading  and  velocity  decay  seen  in 
Figs.  6.8  and  6.9. 

6.5  Conclusions 


In  this  chapter,  axisymmetric  jets  and  wakes  were  studied.  The 
turbulence  model  for  these  flows  is  closely  connected  with  the  plane 
flow  model.  Vortex  stretching  was  identified  as  the  principal  differ¬ 
ence  between  axisymmetric  and  plane  flows.  It  reduces  the  shear  stress 
and,  thereby,  the  integral  parameters  of  the  flow.  Incorporation  of 
this  effect  into  the  model  required  a  change  of  one  model  constant.  The 
readjustment  of  this  constant  was  accomplished  as  In  the  plane  case. 
Excellent  results  for  axisymmetric  co-flowing  jets  and  wakes  were  ob¬ 
tained  with  the  present  model. 

This  model  reduces  to  plane  flow  model  when  there  is  no  vortex 
stretching  (S  -*•  °°)  and  can  be  applied  to  general  three-dimensional 


flows.  Tests  of  this  model  for  three-dimensional  flows  is  therefore 


3 


required.  However,  in  order  to  do  so,  initial  conditions  on  a  plane 
normal  to  the  flow  direction  must  be  provided.  This  requires  a  set  of 
experimental  data  that  covers  the  initial  plane.  Unfortunately,  no 
experimental  data  presently  available  give  this  information.  Most  of 
the  data  were  taken  only  on  two  axes.  For  this  reason,  this  model 


-  *  k  *  -  *  -  •  * 

1  .  ' "!  »  •  , 

» %  ■ y* » » 


cannot  be  tested  for  three-dimensional  flows  at  present.  Measurements 
providing  the  necessary  data  for  Jets  and  for  wakes  of  noncircular  and 
nonplanar  cross  sections  are  accordingly  a  next  step  toward  further 
generalization  of  zonal  nodels  of  free  shear  flows. 


Chapter  VII 


DISCUSSION,  CONCLUSIONS,  AND  RECOMMENDATIONS 

7.1  Summary  and  Discussion 
.  -  .  - *. - 

In  this  research,  we  developed  zonal  models  for  homogeneous  and 
free  shear  flows.  In  each  case,  a  linear  first-order  ordinary  differ¬ 
ential  equation  provides  a  good  readjustment  between  zones.  We  have, 
therefore,  not  explored  more  complex  readjustment  models. 

The  homogeneous  flows  were  modeled  primarily  to  provide  a  quick 
test  of  the  concept  of  zonal  modeling.  Details  of  the  zonal  models  for 
three  types  of  homogeneous  flows  are  given  in  Table  2.1.  The  homogen¬ 
eous  flows  show  that  accurate  prediction  can  be  achieved  by  adjusting 
the  constants  in  the  standard  K-e  model.  The  accuracies  of  the  two 
types  of  modeling  are  indicated  in  Table  7.1;  the  differences  are  based 
on  the  worst  point  on  each  curve.  These  results  partially  validate  the 
OPINION  by  Kline  (1981). 


Table  7.1 


Acccuracy  Comparisons  for  the  Predictions  of  Homogeneous  Flows 


Worst  Point  on  Curve 

Percentage 

of  Cases  Accurate  to 

the  Limit 

Agrees  with  Data  Within 

<  10% 

<  25% 

<  50% 

>  50% 

Zonal  Models 

86 

10 

4 

0 

Standard  K-e  Models 

16 

16 

17 

51 

The  zonal  model  predicts  all  free  shear  flows  within  the  uncer¬ 
tainty  in  the  data.  This  model  is  given  in  detail  in  Table  7.2.  Com¬ 
parisons  of  the  accuracy  of  the  model  of  Table  7.2  with  that  of  the 
standard  K-e  model  are  shown  in  Table  7.3. 
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Table  7.2 

Summary  of  Unified  Zonal  Model  for  Free  Shear  Flows 


Present  Model 

Standard  K-e  Model 

Cu 

0.05  +  0.04  X2  ,  S  2  11-3 

0.09 

c‘l 

1.44  -  0.4  Xj 

1.44 

°e2 

1.92 

1.92 

°K 

1.0  -  0.35  Xj 

1.0 

a 

e 

1.3 

1.3 

Xj  =  1  -  e  CS  ,  X2 


C  =  0.18  ,  C, 


1  -  e 


-C  ( S— 11.3) 
s 


In  general, 


(D61/2/Dt) 

6  '  '  s 
For  steady  flows,  0  simplifies  to: 


si/2 


Table  7.3 

Accuracy  Comparisons  for  the  Predictions  of  Free  Shear  Flows 


Percentage  of  Cases  Accurate  to  the  Limit 


>  10%  >  25%  >  50%  >  50% 


Worst  Point 

on 

Curve 

Agrees  with 

Data  Within 

Unified  Zonal 

Model 

Standard  K 

-c 

Model 

IV  ,F|»T 


As  Tables  7.2  and  7.3  show,  the  zonal  modeling  concept  not  only 
succeeds  but  also  provides  significant  and  unexpected  results. 

First,  while  the  models  apply  to  the  "zones"  of  free  shear  flows, 
it  is  possible  to  connect  these  models  to  yield  a  single  model  contain¬ 
ing  two  non-dimensional  parameters.  This  unified  model  predicts  all  the 
available  free  shear  flow  data,  including  plane  and  axisymmetric  jets 
and  wakes  and  mixing  layers. 

Second,  the  work  enlightened  us  significantly  about  the  underlying 
physics  of  free  shear  flows  and  the  governing  parameters  needed  to 
represent  the  effects  observed. 

These  gains  leave  two  questions  that  need  discussion.  The  seman¬ 
tics  introduced  initially  does  not  describe  the  results  adequately  and 
needs  improvement.  We  also  need  to  know  whether  the  physical  insights 
are  particular  to  this  set  of  flows  or  more  general. 

The  semantic  question  centers  on  the  definition  of  the  word  "zone." 
We  defined  a  zone  to  be  one  or  more  flow  regions  which  could  be  modeled 
by  a  single  model  including  particular  values  of  the  parameters.  In  the 
homogeneous  flows,  this  description  fits  the  situation  precisely;  three 
types  of  strains  needed  three  different  zonal  models.  For  the  free 
shear  flows,  the  situation  is  more  complex.  Zones  representing  the  near 
and  far  fields  were  created  and  linked  using  a  simple  readjustment 
model.  From  that  point  onward,  the  modeling  process  almost  took  on  a 
life  of  its  own.  The  readjustment  from  the  near  to  the  far  field 
employs  a  velocity-ratio  parameter  8.  However,  in  both  pure  jets  and 
mixing  layers  there  is  only  a  single  value  of  this  parameter  (8) 
throughout  the  flow.  We  found  further  that  the  far  field  of  the  axi¬ 
symmetric  jet  is  well  modeled  by  the  plane  far-field  model.  This  made 
it  natural  to  determine  why  the  near  field  of  the  axisymmetric  jet  is 
not  well  modeled.  Vortex  stretching  was  identified  as  the  probable 
cause.  Once  a  parameter,  S,  accounting  for  the  effect  of  vortex 
stretching  was  introduced,  the  unified  model  emerged. 

The  free  shear  flows  do  not  constitute  a  zone  in  the  sense 
Initially  defined.  We  ended  with  three  basic  models  connected  by  two 
parameters.  The  basic  models  (the  physical  situations  they  represent 
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might  best  be  called  "flow  states")  do  not  represent  physical  regions 
but  parts  of  flows  with  limiting  values  of  the  governing  parameters  0 
and  S.  In  flows  in  which  the  parameters  vary,  the  states  are  connected 


by  what  we  originally  called  readjustments  but  might  better  be  called 
"bridges." 


These  ideas  suggest  use  of  two  sets  of  words.  At  an  elemental 


level,  we  might  speak  of  a  zone  and  a  zonal  model  for  that  zone.  Such  an 


elemental  zonal  model  might  still  have  a  readjustment  region  as  from  an 


initial  to  a  final  or  from  a  near  field  to  a  far  field.  At  a  second 


level  of  organization,  such  as  that  of  the  final  results  for  the  class 


of  free  shear  flows,  where  a  single  model  containing  parameters  is 
found,  we  might  speak  of  a  "unified  zonal  model"  for  a  defined  "class". 


What  about  the  question  of  generality  of  the  process?  Inspection 
of  the  processes  used  herein  suggests  that  they  will  generalize  to  at 
least  some  degree.  The  model-construction  process  is  based  on  a  search 
for  inadequacies  in  existing  models,  identification  of  probable  causes 
and  finding  simple  modifications  that  model  them.  This  process  is 
likely  to  succeed  as  well  for  other  flows;  there  is  nothing  special 
about  the  free  shear  flows.  The  latter  were  a  good  choice  for  this 
study,  since  there  is  a  variety  of  cases  in  this  class  of  flows. 


This  augurs  well  for  the  utility  of  zonal  modeling  as  a  method  for 
creating  accurate  models  of  turbulent  flows  for  engineering  purposes.  It 
also  provides  a  tool  for  increasing  our  basic  understanding  of  tech¬ 
nically  significant  turbulent  flows.  Indeed,  the  most  important  result 
of  this  work  may  be  the  process  suggested  rather  than  the  particular 


models  developed. 


7.2  Conclusions 


We  have  shown  that  the  zonal  modeling  concept  can  be  applied  to 
homogeneous  turbulence  and  free  shear  flows.  This  is  a  proof  of  the 
concept  and  provides  an  indication  of  what  might  be  expected  from  zonal 
modeling.  However,  this  is  only  a  very  small  portion  of  a  much  larger 
project  that  needs  to  be  carried  out;  a  great  deal  remains  to  be  done. 
The  zonal  concept  provides  a  logical  framework  within  which  simple  and 
accurate  turbulence  models  can  be  developed  more  easily  than  the 
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traditional  way.  The  present  research  has  laid  a  good  groundwork  for 
the  project  and  we  hope  that  the  encouraging  results  obtained  reflect 
what  can  be  achieved  with  more  complex  flows. 

Finally,  although  we  cannot  anticipate  how  far  this  project  will 
carry  us  toward  reliable  and  fast  computations  of  practically  important 
complex  turbulent  flows,  we  are  confident  that  we  will  be  able  to 
achieve  significantly  more  than  has  been  done  in  the  past. 

7 . 3  Recommendations  for  Future  Work 

The  followings  are  suggestions  for  future  work  in  further  testing 
the  present  unified  model  for  other  free  shear  flows,  understanding  the 
physics,  and  extending  zonal  modeling  to  more  complex  flows. 

1.  A  three-dimensional  free  shear  flow  can  provide  a  challenging  test 
case  for  the  present  unified  model.  However,  detailed  measurements 
on  an  entire  plane  must  be  available  to  start  the  computation.  No 
candidate  flow  has  been  measured  with  sufficient  detail  to  meet 
these  requirements.  Therefore,  experiments  such  as  elliptic  or 
rectangular  jets  with  full  measurements  on  a  cross  section  for  a 
minimum  of  two  stations  would  be  useful. 

2.  The  two  physical  phenomena  parameterized  by  0  and  S  in  the 
present  model  are  mainly  observed  from  the  experimental  results; 
the  causes  of  the  underlying  effects  are  not  understood.  More 
fundamental  research  work,  either  experimental  or  analytical,  that 
can  address  these  causes  might  increase  understanding  of  the  nature 
of  turbulence. 

3.  An  extension  of  the  present  work  to  the  study  of  free  shear  flows 
with  density  difference  and/or  scalar  transport,  such  as  heat  or 
species,  would  make  the  zonal  model  more  valuable  in  practical  en¬ 
gineering  application. 

4.  The  readjustment  parameters  are  solutions  of  first-order  ordinary 
differential  equations  in  the  present  study.  In  later  development 
of  the  zonal  modeling  project,  a  similar  format  in  readjustments 
can  probably  be  used  in  a  complex  flow-field  where  patching  of 
zonal  models  is  required.  This  will  make  the  task  of  creating 
zonal  models  for  complex  flow  fields  easier  and  more  systematic. 
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Fig.  1.1.  Flow  zones  In  a  diffuser  flow  with  separation. 
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Homogeneous  shear  flow  case  2,  dU/dy  “  20  /2,  v  ■  0.02//2. 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo' 
simulation . 
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Homogeneous  axl symmetric  strain  case  1,  dU/dx  =  10,  v  ■  0.01//2 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo's 
slmulat Ion . 


Homogeneous  axlsymmetrlc  strain  case  2,  dU/dx  ”  20,  v  ■  0.01//2 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo's 
simulation . 


Homogeneous  axisyrametric  strain  case  3,  dU/dx  ■  40,  v  *  0.01//2. 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo's 
simulation . 


Homogeneous  axisymmetric  strain  case  5,  dU/dx  *  10,  v  ■  0.02//2. 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo's 
simulation . 


Homogeneous  axl symmetric  strain  case  6,  dU/dx  “  20,  v  ”  0.02//2 
Comparison  of  zonal  model  and  standard  K-e  model  with  Rogallo's 
simulation . 
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Homogeneous  plaice  strain  case  2,  dW/dz  3  1.3, 
n  =  4.299  x  10  Comparison  of  zonal  aodel  and  stan 
dard  K-e  model  with  Lee  and  Reynolds'  simulation. 
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Fig.  2.19.  Homogeneous  pla^e  strain  case  5,  dW/dz  =  25, 

v  ■  4.299  x  10  .  Comparison  of  zonal  model  and  stan 

dard  K-c  model  with  Lee  and  Reynolds'  simulation. 


Homogeneous  plane  strain  case  6,  dW/dz  ■  100, 
v  =  4.299  x  10-  .  Comparison  of  zonal  model  and  stan¬ 
dard  K-e  model  with  Lee  and  Reynolds'  simulation. 


EXPERIMENTAL  DATA 
ZONAL  MOOEL 


STANDARD  K-E  MOOEL 


.22.  Homogeneous  plane  strain,  case  374B  in  80-81  Stanford  Conference, 
dW/dz  »  4.45.  Comparison  of  zonal  model  and  standard  K-e  model 
with  experimental  data. 
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Appendix  A 

FORMULATION  OF  FINITE  DIFFERENCE  EQUATIONS 


Recall  from  Chapter  IV  that  1-0  for  plane  flows  and  i  -  1  for 
axlsyometrlc  flows. 

We  now  consider  the  grid  box  shown  In  Fig.  4.1(b).  We  denote  grid 
points  and  spacings  by 

-  starting  x-location,  £n+*  -  5°  +  A£n,  n  -  0,  1,  2,  ...  (A. 8) 


0;  n. 


n.  +  An.;  j  -  0,  l,  ...,  J;  nT  =  n  (A. 9) 
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and  approximate  the  quantities  (U,V,K, e,T,Q,W)  at  point  (£n,nj)  by 
(Uj,  Vj,  Kj ,  T°,  Qj,  Wj).  We  also  define  the  following  midway 

quantities: 
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(hj  +  nj+1)  (A.  10) 


We  approximate  Gqs.  (A.  1)  to  (A.  4)  at  the  center  of  each  grid  box 
and  Eqs.  (A. 5)  to  (A. 7)  at  midpoint  (£n+*»hj+i/2)*  Keeping  all  the 
known  quantities  (solutions  at  £-  £n)  on  the  right-hand  sides,  we  have 
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In  the  above  equations, 
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Expressions  similar  to  Eq.  (A. 24)  apply  to  0j+i/2»  ^+1/2*  ^j+1/2* 

However,  those  T,  Q,  and  W  associated  with  V  (in  Eqs.  (A.  12)  to 
(A. 14)  have  an  upwind  differencing  switch  in  them  and  are  expressed  as  follows: 
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where  dj+i/2  is  a  switching  parameter  (°j+i/2  *  1  is  ^or  central  differ¬ 
encing  and  Dj+1/2  *0  is  for  upwind  differencing),  and 
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Again,  expressions  similar  to  Eq.  (A.26)  apply  to  Qj+j/2*  Wj+l/2‘  Tj+l/2’ 


O^!  ...  and  Un+!/„  in  those  convection  terms.  The  boundary  conditions  are: 
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